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Granular materials are ubiquitous in our daily lives. While they have been a subject of intensive 
engineering research for centuries, in the last decade granular matter attracted significant attention 
of physicists. Yet despite a major efforts by many groups, the theoretical description of granular 
systems remains largely a plethora of different, often contradicting concepts and approaches. 
Authors give an overview of various theoretical models emerged in the physics of granular matter, 
writh the focus on the onset of collective behavior and pattern formation. Their aim is two-fold: 
to identify general principles common for granular systems and other complex non-equilibrium 
systems, and to elucidate important distinctions between collective behavior in granular and 
continuum pattern-forming systems. 
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iGoUub and Lanee r fl99fj^:',Taege r fiii al\ \l9 9^. 'Kadanof 
"^999); Neddcrman (1992): O ttino and Kh akhar (200qJ; 
iRaichenbacd lj200(l ): iR.istowl l|l999() . Granular materials 
are collections of discrete macroscopic solid grains with 
sizes large enough that Brownian motion is irrelevant (en- 
ergy of 1 mm grain moving with typical velocity of 1 
cm/sec exceeds the thermal energy at least by 10 orders 
of magnitude). Since thermodynamic fluctuations do not 
play a role, for granular systems to remain active they 
have to gain energy either from shear or vibration and are 
thus far from equilibrium. External volume forces (grav- 
ity, electric and magnetic fields) and flows of interstitial 
fluids such as water or air may also be used to activate the 
grains. When subjected to a large enough driving force, a 
granular system may exhibit a transition from a granular 
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solid to a liquid and various ordered patterns of grains 
may develop. Understanding fundamentals of granular 
materials draws upon and gives insights into many fields 
at the frontier of modern physics: plasticity of solids, 
fracture and friction; complex systems from equilibrium 
such as colloids, foams, suspensions, and biological self- 
assembled systems. Moreover, particulate flows are cen- 
tral to a large number of industries including the chemi- 
cal, pharmaceutical, food, metallurgical, agricultural and 
construction industries. Beyond these industrial applica- 
tions, particle laden-flows are widespread in nature, for 
example dune migration, erosion/deposition processes, 
landslides, underwater gravity currents and coastal ge- 
omorphology, etc. 

From a theoretical point of view, it is sometimes use- 
ful to employ an analogy between granular matter and 
ordinary condensed matter and to regard the grains 
as the equivalent of (classical) atoms. However, this 
analogy is far from complete, since the dissipative na- 
ture of grain interactions is the source of many differ- 
ences between the two "kinds" of matter. In partic- 
ular, dissipation is responsible for the fact that most 
states of granular matter are metastable. The typical 
macroscopic size of the grains renders thermal fluctua- 
tions negligible and most standard thermodynamic con- 
cepts inapplicable. Whereas the behavior of dilute gran- 
ular systems (rapid granular gases) can often be ex- 
plained using the f ramework of kinetic theory (see e.g. 
iBrilliantov and P5s chel (2 004j) '). the quantitative theory 
of dense granular assemblies is far less developed. 

In recent years several comprehensive reviews 
and monographs have appea red on the s ubject 
of granular physics, see jAradian all 1200^ ; 
Brilliantov and Poschell20o3:lDuran '. 1999: J aeger etal. , 
19961: iKudrolll 120041: lOttino and Khakhaii l2nnr : 
RaichenbachL l2000t iRistowl l200l|) . Yet in most of them 
the focus has been on actual phenomena and experiments 
rather than on theoretical concepts and approaches to 
the problems of granular physics. Furthermore, the 
scope of granular physics has become so broad that 
we chose to limit ourselves with reviewing the recent 
progress in a subfield of granular pattern formation 
leaving out many interesting and actively developing 
subjects. We loosely define pattern formation as a 
dynamical process leading to the spontaneous emergence 
of nontrivial spatially non-uniform structure which is 
weakly dependent on initial and boundary conditions. 
According to our working definition, we include in the 
scope of the review the patterns in thin layers of vibrated 
grains (Sec. IV, V), patterns in gravity-driven flows (Sec. 
VI), granular stratification and banding (Sec. VII), 
as well as a multitude of patterns found in granular 
assemblies with complex interactions (Sec. VIII). Before 
delving into details of theoretical modelling of these 
pattern- forming systems, we present a brief overview of 
the relevant experimental findings and main theoretical 
concepts (Sec. II and III). 



B. Fundamental microscopic interactions 

Probably the most fundamental microscopic property 
of granular materials is irreversible energy dissipation in 
the course of interaction (collision) between the parti- 
cles. For the case of so-called dry granular materials, 
i.e. when the interaction with interstitial fluid such as 
air or water is negligible, the encounter between grains 
results in dissipation of energy while total mechanical 
momentum is conserved. In contrast to the interaction 
of particles in molecular gases, the collisions of macro- 
scopic grains is generally inelastic. There are several 
well-accepted models addressing the specifics of energy 
dissipation in the course of collision, see for details e.g. 
IBrilliantov and Poschell ("2004). The simplest case corre- 
sponds to in-deformable (hard) frictionless particles with 
fixed restitution coefficient < e < 1 characterizing the 
fraction of energy lost in the course of collision. The 
relation between the velocities after the collisions {v'l 2) 
and before the collision (vi.2) for two identical spherical 
particles is given by 

1 + e 

vi,2 = Vl,2 T [ni2(vi - V2)]ni2. (1) 

Here ni2 is the unit vector pointed from the center of 
particle 1 to the center of particle 2 at the moment of 
collision. The case of e = 1 corresponds to the elastic 
collisions (particles exchange their velocities) and e = 
characterizes fully inelastic collisions. For < e < 1 the 
total energy loss is of the form 

Ai? = -i^|ni2(vi-V2)p. 

Modelling collisions between particles by a fixed restitu- 
tion coefficient is very simple and intuitive, however this 
approximation can be questionable in certain cases. For 
example, approximation of granular media by a gas of 
hard particles with fixed e often yields non-physical be- 
havio r such as inelastic collapse l|McNamara and Yo^md . 
119961) : divergence of the number of collisions in a finite 
time, see Subsec. IIV.AI In fact, the restitution coefficient 
is known to depend on the relative velocities of colliding 
particles and approaches unity as |vi — V2I 0. This 
dependence is captured by the visco-elasti c mode lling of 
particle coUision (see e.g. iRamirez et al\ l)l999|) '). For 
non-spherical grains the restitu tion coefficient m ay also 
depend on the point of contact (Gold smithLll964) . 

Tangential friction forces play an important role in the 
dynamics of granular matter, especially in dense systems. 
Friction forces are hysteretic and history dependent (the 
contact between two grains can be either stuck due to 
dry friction or sliding depending on the history of in- 
teraction). This strongly nonlinear behavior makes the 
analysis of frictional granular materials extremely diffi- 
cult. In the majority of theoretical studies, the simplest 
Coulomb law is adopted: friction is independent on slid- 
ing velocity as long as tange ntial force exceeds the cer- 
tain threshold llWaltonLlT99l . However, the main prob- 
lem is represented by calculation of the static friction 
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forces. It is well known that frictional contact forces 
among solid particles exhibit indeterminacy in case of 
multiple contacts per particle because there are less force 
balance constr aints tha n stress cornponen ts (see, e.g. 
lMcNamara~ g/., (,20o4 . IUnger et all (|2005|) V To resolve 
this indeterminacy in simulations, various approximate 
algorithms have been proposed. In soft particle molecu- 
lar dynamics simulations the most widely used approach 
to calculating friction force s is the spring-dashpot model 
(ICundall and Strackl Il979t ISchaffer et all 1199(1) . An- 
other approach is taken in the contact dynamics method. 
By assuming that all particles are rigid and treating all 
contacting particles as performing instantaneous colli- 
sions (even those which are in fact in persistent contact), 
one can compute the contact forces generated during 
these collisions based on local force balan ce and impen- 
etrabi li ty of the parti cles constraints (see iBrendel et al\ 
l20n4ft : lMoreallll994l) 'l. 

Viscous drag forces due to interaction with intersti- 
tial fluid often affect the dynamics of granular materials. 
Gas-driven particulate flows is an active research area 
in the engineering community, see e.g. Jackson ( 2000). 
Fluid-particle interactions are also invol ved in many geo - 
physical processes, e.g. dune formation l)Bagnoldlll954|) . 
Whereas interaction of small individual particles with 
the fluid is well-understood in terms of Stokes law, col- 
lective interaction and mechanical momentum transfer 
from particles to fluid remains an open problem. Various 
phenomenological constitutive equations are used in the 
enginee ring community to model fluid-particulate flows, 
see e.g. iDuru et all l)2002|) . 

Finally, small particles can acquire electric charge of 
magnetic moment. In this situation fascinating collec- 
tive behavior emerge due to competition between short- 
range colli sions and lo ng-range el ectromagnet i c forces , 
see e.g. lAranson et a l. (200j); iBlair et all l|2003a(l : 
ISapozhnikov et all l)2003j) . Effects of complex inter- 
particle interactions on pattern formation in granular sys- 
tems will be discussed in Sec. I Villi 



II. OVERVIEW OF DYNAMIC BEHAVIOR IN 
GRANULAR MATTER 

In this Section we give a brief overview of the main 
experiments illustrating the dynamical behavior of gran- 
ular media and the phenomena to be discussed in greater 
depth in the following Sections. We classify the exper- 
iments according to the way energy is injected into the 
system: vibration, gravity, or shear. 



FIG. 1 Top view of den se immobile cluster coexist inK with 
dilute granular gas, from lOlafsen and UrbachI (Il998f) . 



FIG. 2 A typical clustering configuration in two dimensions, 
restitution coefficient 0.6, number of particles 40,000, from 
■Goldhirsch and Zanctti (J^Q^J- 



particles, l|01afsen and UrbachL [19981) . Fig. [2 This clus- 
tering transition occurs when the magnitude of vibra- 
tion is reduced (the system is "cooled down") which is 
reminiscent to the clustering instability observed in non- 
driven (freely cooling) gas of inelasti c particles discov- 
ered bv lGoldhirsch and Zanettil l)l993|) . Fig. Detailed 
consideration of clustering phenomena in sub-monolayer 
systems is given in Sec. IIVI 

Multilayers of granular materials subject to vertical vi- 
bration exhibit spectacular pattern formation. In a typ- 
ical experimental realization a layer of granular material 
about 10-30 particle diameters thick is energized by pre- 
cise vertical vibration produced by an electromagnetic 
shaker. Depending on experimental conditions, plethora 
of patterns can be observed, from stripes and squares to 
hexagons and interfaces, see Fig. |21 While the first ob- 
servations of patterns in vi brated layers we re m ade more 
than two centuries ago bv lChladni I l)l787|) and IFaradavl 
(1831,), t he current inte r est in these problems wa s ini- 
tiated by iDouadv et all ljl98£ ):lFauve et all lll989l) and 
culminated in the discovery bv lUmbanhowar et 0,11 lll99fiD 
of a remarkable localized object, oscillon. Fig. ^ De- 
tailed consideration of these observations and their mod- 
elling efforts is given in Sec. |3 

In another set of experiments pattern formation was 
studied in a horizo n tally v i brated system, see e.g. 
iLiffman et all l)l997(l : iRistowl |)1997|) : iTennakoon et all 
(|1998() . While there are certain common features, such 
as sub-harmonic regimes and instabilities, horizontally 
vibrated systems do not show richness of behavior typ- 
ical for the vertically vibrated systems, and nontrivial 
flow regimes are typically localized near the walls. When 
the granular matter is polydisperse, vertical or horizontal 
shaking often leads to segregation. The most well-known 
manifestation of this segregation is the so-called "Brazil 
nut" effect when large particles float to the surface of 
a gra nular layer under vertical shaking ijRosato et all . 
I1987D . Horizontal shaking is also known to produce inter- 
esting segregation band patterns oriented or thogonally to 
the direction of shaking l|MuninLl2000ll2002() (see Fig. 0). 



A. Pattern formation in vibrated layers 



Quasi-two-dimensional sub-monolayers of grains sub- 
jected to vertical vibration exhibit a surprizing bimodal 
regime characterized by a dense cluster of closely packed 
almost immobile grains surrounded by gas of agitated 



FIG. 3 Representative patterns in vertically vibrated granu- 
lar layers for various values of frequency and amplitude of the 
vibration: stripes, squares, hexagons, spiral, interfaces, and 
localized oscillons, from fUmbanhowar et al.. LIQQS) . 
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FIG. 4 Localized oscillon in ver tically vibrate granular layer, 
from iUmbanhowar et ttZ.lllQQ^ . 



FIG. 5 Sequence of snapshots of a layer of copper balls/poppy 
seeds mixture in a horizontally shaken cavity (frequency 12.5 
Hz, amplitude 2 mm) at tim es 5 min, 10 min, 15 min, 30 min, 
1 h, 6 h, from .MullinI (|2000D . 



B. Gravity-driven granular flows 

Gravity-driven systems such as chute flows and 
sandpiles often exhibit nontrivial patterns and spatio- 
temporal structures. Possibly the most spectacular are 
avalanches observed in the layers of granular matter if 
the inchnation exceeds the critical angle (static angle 
of repose). Avalanches were a subject of continued re- 
search for many decades, however only recently it was 
established that the avalanche shape depends sensitively 
on the thickness of the layer and the inclination an- 
gle: triangular downhill avalanches in thin layers and 
balloon-shaped avalanches in thicker layers whi ch ex- 
pand b oth up hill and downhill, see Fig. 1^ and IPaerd 
l)200ll) : Ibaerr and Douadv (199^. Gravity-driven gran- 
ular flows are prone to a variety of non-trivial sec- 
ondary instabilit i es in granular chute flow: fingering 
l|Pouliaiien el ol\ . ^^^), see Fig. E longitudiri a l vor- 
tices in r apid chute flows f'Borzsonvi and Eckel 120051 
iForterre'and Pouhquen. .2001.). see Fig. |^ long modu- 
lation waves (jPorterre and PouliaueiiLl200^ . and others. 

Rich variety of patterns and instabilities has also been 
found in underwater flows of granular matter: transverse 
instability of an avalanche fronts, fingering, pattern for- 
mation in th e sediment behind the avalanche, etc. (see 
iDaerr et al] {2003); Malloeei et al. (2005a b)). Whereas 
certain pattern forming mechanisms are specific to the 
water-granulate interaction, one also finds striking simi- 
larities with the behavior of "dry" granular matter. 



C. Flows in rotating cylinders 

Energy is often supplied into a granular system 
through the shear which is driven by the moving walls 
of the container. One of the most commonly used ge- 
ometries for this class of systems is a horizontal cylin- 
der rotated around its axis, or rotating drum. Rotating 
drums partly filled with granular matter are often used in 
chemical engineering for mixing and separation of parti- 
cles. Flows in rotating drums recently became a subject 



FIG. 6 Sequence of images illustrating evolution of avalanches 
in thin layers on incline. Three images on left: triangular 
avalanche in thin layer, point b in Fig. 1241 Three right im- 
ages: up-hill avalanche in thicker layer, point c in Fig. 1241 
from lPacrr and Douadv Q999il 



FIG. 7 Fingering instability in chute flow, (a) Schematics 
of the instability mechanism, the arrows represent trajectory 
of coarse particles. Images taken from front (b) and bottom 
(c) illustrating accumu lation of coarse particle s between the 
advancing fingers, from IPouliauen et aZ.IJl997D . 



FIG. 8 Development of longitudinal vortices in 
the rapid granular flow down rough incline, from 
iForterre and PouUauenI i200lfl . 

of active research in the physics community. For not too 
high rotating rates the flow regime in the drum is sepa- 
rated into an almost solid-body rotation in the bulk of 
the drum and a localized fluidized layer near the free 
surface (Fig. |5J). Slowly rotating drums exhibit oscilla- 
tions related to the gradual increase of free surface angle 
to the static angle of repose and subsequent fast relax- 
ation to a lower dynamic repose angle via an avalanche. 
Transition t o steady flow is ob served for the higher ro- 
tation rate l|RaichenbachLll990|) . Scaling of various flow 
parameters with the rotation speed (e.g. the width of 
the fluidized layer etc) and development of correlations 
in "dry" and "wet" granular matter was recently studied 
bv lTegzes etai\ ll2002l I200I . 

Rotating drums are typically used to study size 
segregation in binary mixtures of granular materials. 
Two types of size segregation can be distinguished: 
radial and axial. Radial segregation is a relatively 
fast process and occurs after a few revolutions of the 
drum. As a result of radial segregation larger par- 
ticles are expelled to the periphery and a core of 
small er particles is formed in the bulk (Khakhar et ai, 
liggillMetcalfe et a/.LIl995tlMetcalfe and Shattuck.l99£ : 
lOttino and KhakhaiiT200q) . see Fig. HOI 

Axial segregation, occurring in the long drums, hap- 
pens on a much longer time scale (hundreds of rev- 
olutions). As a result of axial segregation, bands of 
segregated materials are formed al ong t he dr um axis 
dHill and Kakahos, 1994, 1995; Zik e Fan . Il994|) . see for 
illustration Fig. ^2 The segregated bands exhibit slow 
coarsening behavior. Even more surprisingly, under cer- 
tain conditions axial segregation pat terns show o scilla- 
tory behavior and trave lling waves ijChoo et alV Il997t 
iFiodor and Ottind. I200I . Possible mechanisms leading 
to axial segregation are discussed in Sec. IVIII 

D. Grains with complex interactions 

Novel collective behaviors emerge when the interac- 
tions between the grains have additional features caused 
by shape anisotropy, interstitial fluid, magnetization or 



FIG. 9 Schematics of flow structure in the cross-section of 
rotating drum, from iKhakhar et al.. (.1997.1 
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FIG. 10 Radial size segregation in a rotating drum, courtesy 
of Wolfgang Losert. 



FIG. 11 Long rotating drum showing axial size segregation, 
from http:/ /www.physics.utoronto.ca/nonlinear/^ 



electrical charge, etc. In this situation short-range col- 
lisions, the hallmark of "traditional" granular systems, 
can be augmented by long-range forces. 

Remarkable patterns including multiple rotating vor- 
tices of nearly ve rtical rods are obser ved in the system of 
vibrated rods bv iBlair et all l)2003a|) . see Fig. [El The 
rods jump on their ends slightly tilted and drift in the 
direction of the tilt. 

Mechanically (" Blair and KudroUiL l20n8H) or electro- 
statically ( Snczh ko et al driven magnetic grains 
exhibit formation of long chains, isolated rings or inter- 
connecting networks, see Fig. ^| In this situation mag- 
netic dipole-dipole interaction augments hard-core colli- 
sions. 

Ordered cl usters and nontrivial d ynami c states were 
observ ed by iThomas and Gollub I l)2004D : IVoth et al\ 
(2002) in a small system of particles vibrated in liq- 
uid (Fig. I14|l . It was shown that fluid-mediated in- 
teraction between particles in a vibrating cavity leads 
to both long-range attraction and short-range repul- 
sion. A plethora of nontrivial patterns including rotat- 
ing vortices, pulsating rings, chains, he xagons etc was 
observed by ISapozhnikov et oil l)2003a|) in the system 
of conducting particles in dc electric field immersed in 
poor electrolyte (Fig. I15|l . The nontrivial competi- 
tion between electrostatic forces and self-induced electro- 
hydrodynamic flows determines the structure of emerging 
pattern. 

Granular systems with complex interactions serve as 
a natural bridge to seemingly different systems such 
as foams, dense colloids, dusty plasmas, ferrofluids and 
many others. 



III. MAIN THEORETICAL CONCEPTS 

Physics of granular media is a diverse and eclectic 
field incorporating many different concepts and ideas, 
from hydrodynamics to the theory of glasses. Conse- 
quently, many different theoretical approaches have been 
proposed to address observed phenomena. 



FIG. 13 Structures formed in submonolayer of magnetic mi- 
croparticles subjected to alternating magnetic field. Select 
structures such as rings, compact clusters, and chains are 
shown in the top panel. Changes in the pattern morphology 
with the increase of magnetic fiel d frequency are illust rated 
by the three bottom images, from lSnezhko et al\ ll2005t) . 

FIG. 14 Regular arrangements of particles near the bottom 
of a vibrated container filled with water when both attraction 
and repulsion are important. All images are taken taken at 
the frequency / = 20 Hz and for different value of dimension- 
less acceleration or for different initial conditions: (a) & (b) 
r = 3; (c) fc (d) r = 3 .7 and F: (e) F = 3.9 and (f) F = 3.5, 
from lVoth etaUi J2002D . 

A. Kinetic theory and hydrodynamics 

Kinetic theory deals with the equations for the proba- 
bility distributions functions describing the state of gran- 
ular gas. The corresponding equations, similar to Boltz- 
mann equations for rarefied gases, can be rigorously de- 
rived for the dilute gas of inelastically colliding parti- 
cles with fixed restitution coefficient, although certain 
generalizations are known , |Gpldstein and Shapiro, 199^ 
iJenkins and ZhanS l2002|) . Kinetic theory is formulated 
in terms of the Boltzmann-Enskog equation for the prob- 
ability distribution function /(v, r, t) to find the particles 
with the velocity v at point r at time t. In the simplest 
case of identical frictionless spherical particles of radius 
d with fixed restitution coefficient e it assumes the fol- 
lowing form 

(at + (vi.V))/((vi,ri,i)=/[/] (2) 
with the binary collision integral /[/] in the form 

/ = (f J dv2 y dni29(-vi2 • ni2)|vi2 • ni2| X 

[x/(vi",ri,t)/(v2",ri-dni2,0 
-/(vi,ri,i)/(v2,ri+dni2,i)] (3) 

where x = l/^^, 6 is theta- function, and pre-coUision 
velocities vi,2 and "inverse collision" velocities z)"2 are 
related as follows 

1 + e 

Vi2 = vi.2T-;^ — [ni2(vi - V2)]ni2 (4) 
2e 

(cf. Eq. (^). This equation is derived with the usual 
"molecular chaos" approximation which implies that all 
correlations between colliding particles are neglected. 



FIG. 12 Select patterns observed in the system of vertically 
vibrated rods with the increase of vibration amplitude: a) 
nematic-like gas phase; b) moving domains of nearly verti- 
cal rod s; c) multiple rotating vortices; d) single vortex, from 
iBlair et al... (.2003a'l . 



FIG. 15 Representative patterns obtained for different values 
of applied field and concentration of ethanol in electrostati- 
cally driven granular system: static clusters (a) and honey- 
combs (b) and dynamic vortices (c) and pulsating rings (d), 
from iSapozhnikov et al.. L2003a) 
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One should keep in mind, however, that in dense granu- 
lar systems this approximation can be rather poor due to 
excluded volume effects and inelasticity of collisions in- 
troducing vel ocity correlati ons among particles (see, for 
example, .Brilliantov and Poschcl (2004,) ). 

Hydrodynamic equations are obtained by truncating 
the hierarchy of moment equations obtained from the 
Boltzmann equation |(2Jl via an appr opriately modified 
Chapman- Enskog procedure (see, e.g.. llBrev ajJ.IT998t 



of state 



iGarzo and Duftvl llQQflt I.Tenkins and Richmanl . ^1985^V 

As a result, a set of continuity equations for mass, 
momentum and fluctuation kinetic energy (or "granu- 
lar temperature" ) is obtained. However, in contrast to 
conventional hydrodynamics, the applicability of gran- 
ular hydrodynamics is often questionable because typi- 
cally there is no separation of scale between microscopic 
and m acroscopic motions ^, see e.g. iTan and GoldhirschI 

The mass, momentum and energy conservation equa- 
tions in granular hydrodynamics have the form 



Dv 

In 

Dm 

'~Dt 
DT 

'~Dt 



-v\I ■ u, 
-V • cr + i/g, 
<j : 7 — V • q — e. 



(5) 
(6) 
(7) 



where v is the filling fraction (the density of granular 
material normalized by the density of grains), u is the 
velocity field, T — ((uu) — (u)^)/2 is the granular tem- 
perature, D/ Dt — dt + {vL-V) is the material derivative, 
g is the gravity acceleration, (Tap is the stress tensor, q is 
the energy flux vector, 7q,/3 = daUp + dpUa is the strain 
rate tensor, and e is the energy dissipation rate. Eqs. 
©-{Tj) are structurally similar to the Navier-Stokes equa- 
tions for conventional fluids except for the last term in 
the equation for granular temperature e which accounts 
for the energy loss due to inelastic collisions. 

These three equations have to be supplemented by the 
constitutive relations for the stress tensor cr, energy flux 
q, and the energy dissipation rate e. For dilute systems, 
a linear relations between stress cr and strain rate 7 is 
obtained, 



<JaP 

q 



[p + (/-t 



A)Tr7](5a/3 - A'7a/3, 



(8) 
(9) 



In the kinetic theory of two-dimensional gas of slightl y 
inelastic hard disks by iJenkins and Richmanl l)l985|) . 
these equations are closed with the following equation 



^ Except the case of almost elastic particles with the restitution 
coefficient e — > 1 



(10) 



and the expressions for the shear and bulk viscosities 



M = 



A = 



l + 2G{v)+ 1 + - G{v) 



7r3/2d ' 
the thermal conductivity 
21.T1/2 



1 + 'iG{v) 



TT^/^dG{y) 
and the energy dissipation rate 

_ l&yG{v)T^''^ 
^- ^3/2^3 



(1 



1 ) G(.) 



(11) 
(12) 

(13) 
(14) 



The radial pair distribution function G{i') for a dilute 
2D gas o f elastic hard disk s can be approximated by the 
formula l|Song eit a/.l . ll989^ 



v{l - 7t//16) 



(15) 



(this is a two-dimensional analog of the famous 
Carn ahan-Starling formula ijCarnahan and StarlineL 
Il969ft for elastic spheres). This formula is expected to 
work for densities roughly below 0.7. For high density 
granular gases, this function has been cal culated using 
free volume theory bv Buehler et al\ l|l95lj) . 



Gpv — 



(16) 



where Vc ~ 0.82 is the d ensity of the random close pack- 
ing limit. iLudind l)200l|) proposed a global fit 



G, 



Gcs + (1 + exp(-(iy - iyo)/mQ))-'^){GFv - Gcs) 



with empirically fitted parameters vq ~ 0-7 and mo ~ 
10^^. However, even with this extension, the contin- 
uum theory comprised of Eas.(Pl- ((ni) cannot describe 
the force chains which transmit stress via persistent con- 
tacts remaining in the dense granular flows, as well as 
the hysteretic transition from solid to static regimes and 
coexisting solid and fluid phases. 

The granular hydrodynamics is probably the most 
universal (however not always the most appropriate) 
tool for modelling large-scale collective behavior in 
driven granular matter. Granular hydrodynamics equa- 
tions in the form ||SJ) , © , Q and their modifications 
are widely used in the engineering community to de- 
scribe a variety of large-scale granular fl ows, espe- 
cially for design of gas-fluidized bed reactors ijGidaspowl 
Il994() . In the physics community granular hydrody- 
namics is used to understand various instabilities in 
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relatively small-scale flows, such as flow past obstacl e 
l|Rericha et a^.l■l2002^■ convectio n llLivn e adl2002alb|) . 
floating chisters llMeerson oil 2003) , longitudinal rolls 
l)Forterre and Pouliaucn'. '2002. 2003), patterns in vi- 
brated layers IBougic et al., 2005) and others. However, 
Eqs. ©-0 are often used far beyond their applicabil- 
ity limits, viz. dilute flows. Consequently, certain pa- 
rameters and constitutive relations need to be adjusted 
heuristically in or der to accommodat e observed behav- 
ior. For example, iBougie et al\ l)2002 lh had to introduce 
artiflcial non-zero viscosity in Eq. ^ for z/ — > in or- 
der to avoid artificial blowup of the solution. Similarly, 
iLosert et al introduced the viscosity diverging as 

density approaches the close packed limit as [u — UcY 
with (3 ~ 1.75 being the fitting parameter in order to 
describe the structure of dense shear granular flows. 



B. Phenomenological models 

A generic approach to t he description of dense granu- 
lar flows was suggested bv lAranson and Tsimrin j l)200ll 
|2p02) who proposed to treat the shear stress mediated 
fluidization of granular matter as a phase transition. For 
this purpose an order parameter characterizing the lo- 
cal state of granular matter and the corresponding phase 
fleld model were introduced. According to the model, the 
order parameter has its own relaxation dynamics and de- 
flnes the static and dynamic contributions to the shear 
stress tensor. This approach is discussed in more details 
in Sec. IVl.A.ll 

Another popular approach is based on the two-phase 
description of granular flow, one phase corresponding to 
rolling grains and the other phase to static ones. This 
approach, so-c alled the BCRE model, was suggested by 
iBouchaud et al. (.1994. .1995.) for description of surface 
gravity driven flows. The BCRE model has direct re- 
lation to depth-averaged hydrodynamic equations (so- 
called Saint- Venant model) popular in the engineering 
community. Note that BCRE and Saint- Venant models 
can be derived in a certain limit from the more general 
order parameter model mentioned above, for detail see 
Sec. IVl.A.21 

Many pattern-forming systems are often described by 
generic amplitude equations such as Ginzburg-Landau or 
Swift-Hohenberg equa t ions llAranson and' Krameill2003 
iCross and Hohenberd. IT991 . This approach allows to 
explain many generic features of patterns, however in 
any particular system there are peculiarities which need 
to be taken into account. This often requires modifl- 
cations to be introduced into the generic model s . Thi s 
approach was taken bv *Aranson and Tsimrinp (1998); 
[Aranson et al. H999a) : Craw ford and Riecke I 199S k 
iTsimring and Aranson I l)l997ij) : IVenkataraniaiiiand Ott I 
l|l99(^ in order to describe patterns in a vibrated gran- 
ular layer. Details of these approaches can be found in 
Sec. rvTDl 

In addition, a variety of tools of statistical physics 



are applied to diverse phenomena occurring in gran- 
ular systems. For example, celebrated theory of 
I Lifshitz and Slvozovl (^58) developed for coarsening 
phenomena in equilibrium systems was successfully ap- 
plied to coarsening of clusters in granular systems 
i)Aranson et"all . l2002() . see Section IVITrni 



C. Molecular dynamics simulations 

Realistic simulation of granular matter consisting of 
thousands of particles remains a challenge for physics and 
computer science. Due to simplicity of microscopic inter- 
action laws (at least for "dry" and non-cohesive granular 
matter) and relatively small number of particles in gran- 
ular flows as compared to atomic and molecular systems, 
the molecular dynamics simulations or discrete element 
models have a potential to address adequately many phe- 
nomena occurring in the granular systems. 

There exist three fundamentally different approaches, 
so-called soft particles simulation method; event driven 
algorithm and the contact dynamics method for rigid 
particles. For the review on various m olecular dynam - 
ics simulation methods we recommend 'Ludind (|200J); 



lind 



iPoschel and Schwager (,2005j) : iRaoaport 1,199£ 

In the soft particle algorithm, all forces acting on a 
particle either from walls or other particles or external 
forces are calculated based on the positions of the par- 
ticles. Once the forces are found, the time is advanced 
by the explicit integration of the corresponding Newton 
equations of motion. Various models are used for cal- 
culating normal and tangential contact forces. In ma- 
jority of implementations, the normal contact forces are 
determined from the particle overlap A„ which is de- 
flned as the difference of the distance between the cen- 
ters of mass of two particles and the sum of their radii. 
The normal force F„ is either proportional to A„ (lin- 
ear Hookian contact) or proportional to A"^/^ (Hertzian 
contact). In the spring-dashpot model, additional dis- 
sipative force proportional to the normal component of 
the relative velocity is added to model inelasticity of 
grains. A variety of approaches are used to model tan- 
gential forces, the most widely accepted of them be ing 
Cundall-Strack algorithm l)Cundall and Stra'cilll979j) . in 
which the tangential contact is modelled by a dissipa- 
tive linear spring whose force Ft = — /cjAj — mjl'^t^t 
(here Aj is the relative tangential displacement and vj 
is the relative tangential velocity, fcj , 7t are model con- 
stants). It is truncated when its ratio to the normal 
force |rt|/|F„| reaches the friction coefficient /i accord- 
ing to the Coloumb law. Soft-particles methods are rela- 
tively slow and used mostly for the analysis of dense flows 
when generally faster event-driven algorithms;^ are not 
applicabl e, see e.g. |Laiid ry_ela^ l)2003j) : ISilbert et al\ 
ll2nn2alblOl : IVolfson et oi.l ll2003albD . 

In the event-driven algorithm, the particles are con- 
sidered inflnitely rigid and move freely (or driven by 
macroscopic external flelds) in the intervals between (in- 
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stantaneous) collisions. The algorithm updates veloci- 
ties and positions of the two particles involved in a bi- 
nary collision (in the simplest frictionless case, accord- 
ing to Eq. (Q), and then finds the time of the next 
collision and velocities and positions of all particles at 
that time according to Newton's law. Thus, the time 
is advanced directly from one collision to the next, and 
so variable time step is dictated by the interval be- 
tween the collisions. While event-driven methods are 
typically faster for dilute rapid granular flows, they be- 
come impractical for dense flows where collisions are 
very frequent and furthermore particles develop persis- 
tent contacts. As a related numerical problem, event- 
driven methods are known to suffer from so-called "in- 
elastic collapse" when the num ber of collisions between 
particles diverges in finite time l)McNamara and Yound . 
[T996). There are certain modifications to this method 
which allow to circumvent this problem by intro- 
duci ng velocity d epende nt restitution cocfhcient (see, 
e.g. iBizon et all l)l998aj) '). but still event-driven meth- 
ods are mostly applied t o rapid granular flows, see 
e.g. iFerguson et a/.l J2004): iKhain an d Mcerson ( 2004,1 : 
iMcNamara and Yound l|l996(l : iNie et ali l|2002D . 



Contact dynamics is a discrete element method like 
soft-particles and event-driven ones, with the equations 
of motion integrated for each particle. Similarly to event- 
driven algorithm and unlike soft-particles method, parti- 
cle deformations are suppressed by considering particles 
infinitely rigid. The contact dynamics method considers 
all contacts occurring within a certain short time interval 
as simultaneous, and computes all contact forces by sat- 
isfying simultaneously all kinematic constraints imposed 
by impenetrability of the particles and the Coulomb fric- 
tion law. Imposing kinematic constraints requires con- 
tact forces (constraint forces) which cannot be calcu- 
lated from the positions and velocities of particles alone. 
The constraint forces are determined in such a way 
that constraint-violating accelerations are compensated. 
For comprehensive r eview on the contact dynamics see 
iBrendel et ai\ l|2004D . 



Sometimes different molecular dyn amics methods ar e 
often applied to th e same problem iLois et al\ ()2005(l : 
iRadiai etHH ()l998(l : IStaron et al\ (l2002t) applied con - 
tact dynamic s method s and ISilbert et ali l|2002albl) : 
IVolfson et al\ ()2003allJ ) used soft-particles technique 
for the analysis of instabilities and constitutive re- 
lations in dense granular systems. Patterns in vi- 
br ated layers w ere studied by event - driven simulations 
bv iBizon et"a?] (|l998a); Moon et al\ l)2003|) a nd by soft 
partic les molecular dynamics simulations by iNie et all 
l|2nnnD : lPrevost fit a.l\ (I2004D . 



IV. PATTERNS IN SUB-MONOLAYERS. CLUSTERING, 
COARSENING AND PHASE TRANSITIONS 

A. Clustering in Freely Cooling Gases 

Properties of granular gases are dramatically different 
from the properties of molecular gases due to inelastic- 
ity of collisions between the grains. This leads to the 
emergence of correlation between colliding particles and 
violation of the molecular chaos approximation. This 
in turn gives rise to various pattern-forming instabili- 
ties. Perhaps the simplest system exhibiting nontrivial 
pattern formation in the context of granular matter is 
freely cooling granular gas: isolated system of inelasti- 
cally colliding particles. The interest to freely cooling 
granular gases was triggered by th e d iscovery of clus- 
tering bv iColdhirsch and Zanettil l)l993j) : spontaneously 
forming dense clusters emerge as a result of instability 
of initially homogeneous cooling state, see Fig. |21 This 
instability, which can be traced in many other granular 
systems, has a very simple physical interpretation: local 
increase of the density of granular gas results in the in- 
crease in the number of collisions, and, therefore, further 
dissipation of energy and decrease in the granular tem- 
perature. Due to proportionality of pressure to the tem- 
perature, the decrease of temperature will consequently 
decrease local pressure, which, in turn, will create a flux 
of particles towards this pressure depression, and further 
increase of the density. This clustering instability has in- 
teresting count erparts in astrophysics: clusterin g of self- 
gravitating gas l|Shandarin and ZeldovichLll989l) and "ra- 
diativ e instability" in optically thin plasmas ijMeersonL 
Il996l) resulting in interstellar dust condensat ion. 

According to lGoldhirsch and Zanettil l)l993() , the initial 
stage of clustering can be understood in terms of the 
instability of a homogeneously cooled state described by 
the density v and granular temperature T. This state is 
characterized by zero hydrodynamic velocity v, and the 
temperature evolution follows from the energy balance 
equation 

dtT - -r3/2 (17) 

which results in the Haff's cooling law T ^ t~'^ (Hafl| 
Il983j) . However, the uniform cooling state becomes 
unstable in large enough systems masking Haff's law. 
The discussion o f the l i near i n stability conditions can 
be found e.g. in iBabid l)l993() : IBrilliantov and Pdschell 
l|2004|) . For the case of particles with fixed restitu- 
tion coefficient e, the analysis in the framework of lin- 
earized hydrodynamics equations ©-((Tl) yields the criti- 
cal wavenumber k* for the clustering instability 

k* - Vl-e2 (18) 

As one sees, the length scale of the clustering instability 
diverges in the limit of elastic particles e ^ 1. 

The clustering instability in a system of grains with 
constant restitution coefficient results in the inelastic 
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collapse discussed in the previous Section. Whereas 
the onset of clustering can be well-understood in 
the framew' o rk of granular hydrodynamic s (see, e.g. 
iBabic' ('iggS"); BriUiantov and Poschcl (2004');"Goldhirsch' 
(|2003. ) : Hill and Mazenko. (^2003.) ). certain subtle features 
(e.g. scaling exponents for temperature) are only as- 
sessed within molecular dynamics simulation because the 
hydrodynamic description often breaks in dense cold clus- 
ters. One recent theoretical approach to the description 
of the late stages of clustering instability consists in intro- 
ducing additional regularization into the h ydrodynamic 
description due to the finite size of particles jEfrati et all 
120051: iNie fi^ oiJ. 120021) . 

iNiefiT^iJ l|2002^ "argued that cluster formation and co- 
alescence in freely cooling granular gases can be heuris- 
tically described by the Burgers equation for hydrody- 
namic velocity v with random initial conditions: 

dtv + vVv = fioV'^v (19) 

where is effective viscosity (which is different from the 
shear viscosity in hydrodynamic equations). In this con- 
text clustering is associated with the formation of shocks 
in the Burgers equation. Perhaps not surprising, a very 
similar approach was applied for description of the gas 
of "sticky" particles for the descri ption of the large-scale 
matter formation in the Univers e ijClurbatov et oi.Lfl985t 
IShandarin and Zeldovich, 1989). 

.Meerson and Puglisi (2,005) conducted molecular dy- 
namics simulations of the clustering instability of a freely 
cooling dilute inelastic gas in a quasi-one-dimensional set- 
ting. This problem was also examined in the fra mework 
of granular hydrodynamics bv lEfratT et a/.ll|2005|) . It was 
observed that, as the gas cools, stresses become negligibly 
small, and the gas flows only by inertia. Hydrodynamic 
description reveals a finite-time singularity, as the veloc- 
ity gradient and the gas density diverge at some location. 
The molecular dynamics studies show that finite-time 
singularities, intrinsic in such flows, are arrested only 
when close-packed clusters are formed. It was confirmed 
that the late-time dynamics and coarsening behavior are 
describable by the Burgers equation (|19|l with vanishing 
viscosity /ig. Correspondingly, the average cluster mass 
grows as t^^^ and the average velocity decreases as t~^/^. 
Due to the clustering long-term temperature evolution is 
T - ^j^jpjj jg different from Haff's law r ^ de- 

rived for the spatially-homogeneous cooling. lEfrati et al\ 
l)2005f) argue that flow by inertia represents a generic in- 
termediate asymptotic of unstable free cooling of dilute 
granular gases consistent with the Burgers equation H19|l 
description o f one-dim e nsion al gas of "sticky particles" 
suggested bv lNie et ai\ l)2002() . 

While there is a qualitative similarity between Burgers 
shocks and clusters in granular materials at least in one 
dimension, the applicability of the Burgers equation for 
the description of granular media is still an open ques- 
tion, especially in two and three dimensions. The main 
problem is that the Burgers equation can be derived from 
the hydrodynamic equations only in one dimensional sit- 



uation, in two and three dimensions the Burgers equation 
assumes zero vorticity, which possibly oversimplifies the 
problem and may miss important physics. In fact, molec- 
ular dynamics simulations illustrate the development of 
large-s cale vortex flows in the cou r se of clustering insta- 
bility l)Catuto and MarconiL 12004 Ivan Noiie and Ernsl 
120001) . 

Finally, iDas and PurJ l)2003tl proposed a phenomeno- 
logical description of the long-term clusters evolution in 
granular gases. Using the analogy between clustering 
in granular gases and phase-ordering dynam ics in two- 
component mixtures, iDas and PuTT |)2003|) postulated 
generalized Cahn-Hilliard equations for the evolution of 
density u an d complex velocity t/j — + ivy (see e.g. 
llB^[l99l ) 

dtu = (-V^)™ [z/-z/3^v2i/] (20) 

dt^ = (-V^)™ [i/.- |V;|V + VV] (21) 

with TO 0+ which characterize globally-conserved dy- 
namics of ly and ^ simi l ar to that considered in Sec. 
IVIII.C.ll iDas and Puril l)2003|) argue that this choice 
is most appropriate due to the non-diffusive character 
of particles motion and is consistent with the observed 
morphology of clusters. While it might be very challeng- 
ing to derive Eqs. (|20|l . (|21|) from the first principles or 
to deduce them from hydrodynamic equations, the con- 
nection to phase-ordering dynamics is certainly deserves 
further investigation. 

B. Patterns in Driven Granular Gases 

Discovery of the clustering instability stimulated a 
large number of experimental and theor etical studies , 
even experiments in low gravity conditions l)Falcon et al\ . 
I1999D . Since "freely cooling granular gas" is difficult to 
implement in the laboratory, most experiments were per- 
formed in the situation when the energy is injected in the 
granu lar system in one or another way. KudroUi et all 
l)l997j) studied two-dimensional granular assemblies in- 
teracting with a horizontally vibrating (or "hot") wall. 
In agreement with granular hydrodynamics, maximum 
gas density occurs opposite to the vibrating wall, see 
Fig. 1161 The experimental density distributions are con- 
sistent wi th the modified hydrod y namic approach pro- 
posed by iGross man et "aT (1997). Khain and M eersonI 
i)2002|) : lKhain et aL (,2004a) : .Livne et aL (2002a.bi) . stud- 
ied the dynamics of granular gases interacting with a hot 
wall analytically using granular hydrodynami c theory for 
rigid disks in the formulation of Jenkins and RichmanI 
l)l985j) and predicted a novel phase-separation or van der 
Waals-type instability of the one-dimensional density dis- 
tribution. This instability , reproduced later by m olecu- 
lar dynamics simulations ijArgentina et all |2002|) is dif- 
ferent from the usual convection instability as it occurs 
without gravity and is driven by the coarsening mecha- 
nism. Simulations indicated a profound role of fluctua- 
tions. One may expect that noise amplification near the 
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FIG. 16 Sample image showing dense cold cluster formed op- 
posite the drivin g wall (at the botto m), total number of par- 
ticles 1860, from lKudrolh et al\E997i) . 



instability thresholds in granular systems wiU be very 
important due to non-macroscopic numb er of grains. In 
the context of phase-separation instabihtv lMeerson et al\ 
12004) raised the non-trivial question of the origin of gi- 
ant fluctuations and break-down of hydrodynamic de- 
scription in granular syste ms ne ar the threshol d of in- 
stabi lity (see also ( Bougie 'ei all |2005; Goldma n et all 
I2004D on the effect of fluctuations in multilayers). Re- 
markably, for the granular gas confined be tween two os- 
cillating walls iKh ain and MeersonI l)2004() predicted on 
the basis of event-driven simulations a novel oscillatory 
instability for the position of the dense cluster. These 
predictions, however, have not yet been confirmed exper- 
imentally, most likely due to the relatively small aspect 
ratio of available experimental cells. 

lOlafsen and UrbachI l)l998|) pioneered experiments 
with sub-monolayers of particles subject to vertical vibra- 
tion^. Their studies revealed a surprising phenomenon: 
formation of a dense closely-packed cluster co-existing 
with dilute granular gas, see Fig. ^ The phenomenon 
bears a strong resemblance to the first-order solid/liquid 
phase transi tion in equil ibrium systems. Similar ex- 
periments bv lLosert et al discovered propagating 
fronts between gas-like and solid-like phases in vertically 
vibrated sub-monolayers. Such fronts are expected in 
extended systems in the vicinity of the first order phase 
transi tion, e.g. solidificat ion fronts in supercooled liq- 
uids. iPrevost et al\ l)2004^ performed experiments with 
vibrated granular gas confined between two plates. Qual- 
itatively similar phase coexistence was found. The cluster 
formation in vibro-fluidized sub-monolayers shares many 
common features with processes in freely cooling gran- 
ular gases because it is also caused by the energy dissi- 
pation due to inelasticity of collisions. However, there is 
a significant difference: the instability described in Sub- 
sec. lIV.Al is insufficient to explain the phase separation. 
A very important additional factor is bistability and co- 
existence of states due to the nontrivial density depen- 
dence of the transfer rate of particle's vertical to horizon- 
tal momentum. Particles in a dense closed-packed cluster 
likely obtain less horizontal momentum than in a moder- 
ately dilute gas because in the former particle vibrations 
are constrained to the vertical plane by interaction with 
neighbors. In turn, in a very dilute gas the vertical to 
horizontal momentum transfer is also inhibited due to 
lack of particle collisions. Another factor here is that vi- 
bration is not fully equivalent to the interaction with a 
heat bath. It is well known that even a single particle 



^ Sub-monolayer implies less than 100% percent coverage by par- 
ticles of the bottom plate. 



interacting with a periodically vibrating plate exhibits 
coexi stence of dynamic and static states ijLosert et all 
I1999D . 

There were several simulation studies of clustering 
and pha se coexist e nce in vibrated gran u lar su bmono- 
layers. iNie et al\ l|2000() : iPrevost et oil l|2004|) repro- 
duced certain features of cluster formation and two-phase 
co-existence by means of large-scale three-dimensional 
molecular dynamics simulations. Since realistic three- 
dimensional simulations are still expensive and extremely 
time-consuming, simplified modelling of the effect of a 
vibrating wall by a certain multiplicative random forc- 
ing on individual particles was employed bv lCafiero et all 
|l99i| ) . While the multiplicative random forcing is an in- 
teresting theoretical idea, it has to be used with caution 
as it is not guaranteed to reproduce subtle details of par- 
ticle dynamics, especially the sensitive dependence of the 
vertical to horizontal momentum transfer as the function 
of the density. 



C. Coarsening of clusters 

One of the most intriguing questions in the con- 
text of phase coexistence in vibrofiuidized granular 
sub-monolayers is a possibility of Ostwald-type ripen- 
ing and coarsening of clusters simila r to that ob- 
served in equilibrium systems ( Lifs hitz and Slvozovl 
Il958l IT96li) . In particular, the scaling law for the num- 
ber of macroscopic clusters is of special interest because 
it gives a deep insight into the similarity between equi- 
librium thermodynamic systems and no n-eq uilibrium 
granu l ar systems. The exper i ments (^ osert et a^. 
19991 lOlafsen and UrbachL Il998t IPrevost et a l. 200^ 
SapozIinikoT^r^in200^r demonstrated emergence and 
growth of multiple clusters but did not have sufficient 
aspect ratio to address the problem of coarsening in a 
quantitative way. 

Nev ertheless, as it was suggested by lAranson et all 
l)2000() , statistical information on out-of-equilibrium Ost- 
wald ripening can be obtained in a different granular sys- 
tem: electrostatically driven granular media. This sys- 
tem permits one to operate with extremely small par- 
ticles and obtain a very large number of macroscopic 
clusters. In this system the number of clusters N de- 
cays with time as ~ 1/t. This law is consistent 
with interface-co ntrolled Ostw ald ripening in two di- 
mensions, see ( W agneij Il96l|) . Whereas mechanisms 
of energy injections are different, both vibrofiuidized 
and electrostatically-driven systems show similar behav- 
ior: macroscopic phase separation, coarsening, transi- 
tion from two- to three-dimensional cluster growth, etc 
fSapozhnikov et at, 2003). In Aranson et al. (2002) the 
theoretical description of granular coarsening was devel- 
oped in application to the electrostatically driven grains, 
however we postpone the description of this theory to 
Sec. IVIII.CI We anticipate that a t heory similar to that 
formulated in lAranson et al ](|200l can be applicable to 
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mechanically fluidized granular materials as well. The 
main difference there is the physical mechanism of en- 
ergy injection which will possibly affect the specific form 
of the conversion rate function (j) in Eq. H()9|l in Sec. 

iviim 



V. SURFACE WAVES AND PATTERNS IN VIBRATED 
MULTILAYERS OF GRANULAR MATERIALS 

A. Chladni patterns and heaping 

Driven granular systems often manifest collective fluid- 
like behavior: shear flows, co nvection, surf a ce wa ves, and 
pattern formation (see e.g. iJaeger et all l)l996j) '). Sur- 
prisingly, even very thin (less than ten) layers of sand 
under excitation exhibit pattern formation which is quite 
similar (however with some important differences) to the 
corresponding patterns in fluids. One of the most fasci- 
nating examples of these collective dynamics is the ap- 
pearance of long-range coherent patterns and localized 
excitations in vertically-vibrated thin granular layers. 

Experimental studies of vibrated layers of sand have a 
long and illustrious histor y, be ginning f r om th e seminal 
works by .Chladni i l|l787|) and iFaradavl l|l83l|) in which 
they used a violin bow and a membrane to excite vertical 
vibrations in a thin layer of grains. The main effect ob- 
served in those early papers, was "heaping" of granular 
matter in mounds near the nodal lines of the membrane 
oscillations. This behavior was immediately (and cor- 
rectly) attributed to the "acoustic streaming" , or nonlin- 
ear detection of the nonuniform excitation of grains by 
membrane modes. One puzzling result by Chladni was 
that a very thin powder would collect at the anti-nodal 
regions where the amplitude of vibrations is maximal. As 
Faraday demonstrated by evacuating the container, this 
phenomenology is caused by the role of air permeating 
the grains in motion. Evidently, the interstitial gas be- 
comes important as the terminal velocity of a free fall 
Vt = fgcP/ 18 fj, becomes of the order of the plate velocity, 
and this condition is fulfilled for 10 — 20 fim particles on 
a plate vibrating with frequency 50 Hz and acceleration 
amplitude g. 

In subsequent years the focus of attention was di- 
verted from dynamical properties of thin layers of vi- 
brated sand, and only in the last third of the 20th cen- 
tury physicists returned to this old problem equipped 
with new experimental capabilities. The dawn of the 
new e ra was marked by the stu dies of heaping by I Jenny! 
lll964l). In subsequent papers llDinkelacker et ad Il987 : 
Douadv et a j. .119891 lEvesaue et a/.l . ll989HLaroche et al . 
1989HWalkeilll982D . more research has been performed of 
heaping with and without interstitial gas, with somewhat 
controversial conclusi ons on the nece ssity of ambient gas 
for heaping (see, e.g. llEvesauel Il990l)') . Eventually, after 
more careful analysis iPak et al\ l)l995|) concluded that 
heaping indeed disappears as the pressure of the am- 
bient gas tends to zero or the particle size increases. 



This ag reed with numerical m o lecular dynamics simu- 
lations (ICaflas et all Il992albll(l | Gallas and Sokolowskil 
119931: iLuding et all 11994 iTaeuchl Il992|) which showed 
no heaping without interstitial gas effects. Recent stud- 
ies of deep layers (50 < N < 200 ) of s r nall particles 
10 < d < 200iim) bv lDuranI l|2000l 12001(1 : iFalcon et all 
1999a^ showed a number of interesting patterns and 
nov el instabilitie s caused by interstitial air. In particu- 
lar, |Dura3 l|200l|) observed formation of isolated droplets 
of grains after periodic taping similar to the Rayleigh- 
T aylor instability in ordinary fiuids. 

Ijia et al\ l)l999|) proposed a simple model for heap for- 
mation which is motivated by these experiments. In a 
discrete lattice version of the model, the decrease in lo- 
cal density due to vibrations is modelled by the random 
creation of empty sites in the bulk. The bulk flow is sim- 
ulated by the dynamics of empty sites, while the surface 
flow is modelled by rules similar to the sandpilc model 
(see Sec. IVI.B|I . This model reproduced both convection 
inside the powder and the heap formation for sufficiently 
large probability of empty sit e formation (whi ch mimics 
the magnitude of vibration). iJia et al] l|l999l) also pro- 
posed the continuum model which has a simple form of a 
nonlinear reaction-diffusion equation, for the local height 
of the sandpile 



dth = D\/^h + nh- (3h^ 



(22) 



However this model is perhaps too generic and lacks the 
specific physics of the heaping process. 



B. Standing wave patterns 

While heaping may or may not appear depending on 
the gas pressure and the particle properties at small ver- 
tical acceleration, at higher vertical acceleration patterns 
of standing waves emerge i n thin layers. They were fi rst 
reported bv 'Douad v et al\ l|l989|) : iFauve et al\ l|l989() in 
a quasi two-dimensional geometry. These waves oscil- 
lated at the half of the driving frequency, which indi- 
cates the sub-harmonic resonance characteristic for para- 
metric instability. This first observation spurred a num- 
ber of experimental studies of standing waves in thin 
granu lar layers in two and t hree dimension al geometries 
(Aran son eroH., 199 9b: Cle ment aLLlT996: Melo et al. 
1994, 199M iMuiicaTand Meld Il998l: lUmbanhowar et al. 



199fiil . Importantly, these studies were performed in evac- 
uated containers, which allowed to obtain reproducible 
results not contaminated by heaping. Fig. |31 shows a 
variety of regular patte rns observed in vi brated granular 
layers under vibration l|Melo et all Il994|) . As a result of 
these studies, the emerging picture of pattern formation 
appears as follows. 

The particular pattern is determined by the interplay 
between driving frequency / and acceleration of the con- 
tainer r = An^Af'^/g {A is the amplitude of oscil l ations , 
g is the gravity acceleration) <|Melo et all 11994 Il995fl . 
The layer of grains remains flat for T < 2.4 more-less 
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FIG. 17 Phase diagram of various regimes in vibrated gran- 
ular layers, from lMelo et al\ ^^9^. 

independent of driving frequency. At higher T patterns 
of standing waves emerge. At small freq uencies f < f* 
(for experimental conditions of lMelo et a l. (1995) , /* ~ 
45 Hz) the transition is subcritical, leading to the for- 
mation of square wave patterns, see FiglJls. For higher 
frequencies / > /* the selected pattern is quasi-one- 
dimensional stripes (Fig. and the transition be- 
comes supercritical. In the intermediate region / ^ /*, 
localized excitations (oscillons, FigQ and various bound 
states of oscillons (FigJSf) were observed within the hys- 
teretic region of the parameter plane. Both squares and 
stripes, as well as oscillons, oscillate at the half of the 
driving frequency, which indicates the parametric mech- 
anism of their excitation. The wavelength of the cellular 
patterns near the onset scales linearly with the depth of 
the laver a nd diminishes with the frequency of vibration 
l|Umbanhowar and Swinnevl l2000l|) . The frequency cor- 
responding to the strip-square transition was shown to 
depend on the particle diameter d as c?^^/^. This scaling 
suggests that the transition is controlled by the relative 
magnitude of the energy influx from the vibrating plate 
oc and the gravitational dilation energy oc gd. At 
higher acceleration (F > 4), stripes and squares become 
unstable, and hexagons appear instead (Fig. |3|;). Further 
increase of acceleration at F « 4.5 converts hexagons into 
a domain-like structure of flat layers oscillating with fre- 
quency //2 with opposite phases. Depending on parame- 
ters, interfaces separating flat domains, are either smooth 
or "decorated" by periodic undulations (Fig. For 
r > 5.7 various quarter-harmonic patterns emerge. The 
complete phase diagram of different regimes observed in 
a three-dimensional container is shown in Fig. 1171 For 
even higher acceleration (F > 7) the experiments reveal 
surprising phase bubbles and spatio-temporal chaos oscil- 
lating approximate ly at one fourth the driving frequency 
l)Moon et a/.Ll2002|) . 

Subsequent investigations revealed that periodic pat- 
terns share many features with convective rolls in 
Rayleigh-Benard convecti on, for example skew-v aricose 
and cross-roll instabilities ijde Bruvn et all Il998() . 



C. Simulations of vibrated granular layers 

The general understanding of the standing wave pat- 
terns in thin granular layers can be gained by the anal- 
ogy with ordinary fluids. The Faraday instability in flu- 
ids and corresponding pattern selection problems have 
been studied theoretically and numeric ally in great de- 
tail (see e.g. IZhang and Vinalj l)l997j) ). The primary 
mechanism of instability is the parametric resonance be- 
tween the spatially uniform periodic driving at frequency 
/ and two counter-propagating gravity waves at fre- 
quency f/2. However, this instability in ordinary flu- 



FIG. 18 Comparison between subharmonic patterns in exper- 
iment (left) and three dimensional molecular dynamics sim- 
ulations (right) of 30000 particles in a square vibrated con- 
taine r for different frequ encies and amplitudes of vibration, 
from lBizon et al\ lll998at) . 



ids leads to a supercritical bifurcation and square wave 
patterns near offset, and as a whole the corresponding 
phase diagram lacks the richness of the granular sys- 
tem. Of course this can be explained by the fact that 
there are many qualitative differences between granu- 
lar matter and fluids, such as presence of strong dis- 
sipation, friction and the absence of surface tension in 
the former. Interestingly, localized oscillon-type objects 
were subsequently observed in vertically vibrated lay- 
ers of non-Newtonian fluid ( Lioubashevski ^J aL, 1999), 
and stipe patterns wer e observed in highly viscous fluid 
l)Kivashko all Il996(l . The theoretical understanding 
of the pattern formation in a vibrated granular system 
presents a challenge, since unlike fluid dynamics there 
is no universal theoretical description of dense granular 
flows analogous to the Navier-Stokes equations. In the 
absence of this common base, theoretical and computa- 
tional efforts in desc ribing these patter ns followed several 
different directions. lAoki et all i)l996(l were first to per- 
form molecular dynamics simulations of patterns in the 
vibrated granular layer. They concluded that grain-grain 
friction is necessary f or pattern formatio n in this system. 
However, as noted bv lBizon et al ] |ll997ft . this conclusion 
is a direct c onseq uence of the fact that the algorithm of 
lAoki et all l)l996j) . which is based on the Lennard- Jones 
interaction potential and velocity dependent dissipation, 
leads to the restitution coefficient of particles approach- 
ing unity for large collision speeds rather than decreasing 
according to experiments. 

iBizon et al\ l)l998albt) performed event-driven simula- 
tions of colliding grains on a vibr ated plat e assuming con- 
stant restitution (see also Luding et al\ jl996 ) for ear- 
lier two-dimensional event-driven simulations). It was 
demonstrated that even without friction, patterns do 
form in the system, however only supercritical bifurca- 
tion to stripes is observed. It turned out that friction 
is necessary to produce other patterns observed in ex- 
periments, such as squares and //4 hexagons. Simula- 
tions with frictional particles reproduced the majority of 
patterns observed in experiments and many features of 
the bifurcation diagram (wi t h the im portant exception of 
the oscillons). iBizon et 

Ml il998a') set out to match an 
experimental cell and a numerical system, maintaining 
exactly the same size container and sizes and the num- 
ber of particles. After fitting only two parameters of the 
numerical model, Bizon et al. ( 1998ai) were able to find 
a very close quantitative agreement between various pat- 
terns in the experimental cell and patterns in simulations 
throughout the parameter space of the experiment (fre- 
quency of driving, amplitude of acceleration, thickness of 
the layer), see Fig. ^1 
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IShinbro^ l)l997() proposed a model which combined 
ideas from molecular dynamics and continuum mod- 
elling. Specifically, the model ignored vertical compo- 
nent of particle motion and assumed that impact with 
the plate adds certain randomizing horizontal velocity to 
the individual particles. The magnitude of the random 
component being added at each impact served as a mea- 
sure of impact strength. After the impact particles were 
allowed to travel freely in the horizontal plane for a cer- 
tain fraction of a period after which they inelastically 
collide with each other (a particle acquires momentum 
averaged over all particles in its neighborhood). This 
model did reproduce a variety of patterns seen in experi- 
ments (stripes, squares, and hexagons) for various values 
of control parameters (frequency of driving and impact 
strength), however it did not describe some of the ex- 
perimental phenomenology (localized objects as well as 
interfaces) , besides it also produced a number of intricate 
patterns not seen in experiments. 

D. Continuum theories 

The first continuum models of pattern formation 
in vibrating sand were purely phenomenological. In 
the spirit of wea kly-nonlinear perturbation theories 
iTsimring and Aranson (.1997.) introduced the complex 
amplitude ip{x,y,t) of sub-harmonic oscillations of the 
layer surface, h = V' exp(z7r/t) -I- c.c. The equation for 
this function on the symmetry grounds in the lowest or- 
der was written as 

dtip = 7V'* - (1 - iuj)^ + (1 + i6)V V - IV'I V - (23) 

Here 7 is the normalized amplitude of forcing at the driv- 
ing frequency /. The linear terms in Eq. 1231) can be 
obtained from the complex growth rate for infinitesimal 
periodic layer perturbations h ~ cxp[A(k)t + ikx]. Ex- 
panding A(fc) for small fc, and keeping only two leading 
terms in the expansion A(fc) = — Aq — Aifc^ gives rise 
to the linear terms in Eq. (|23|l . where b = ImKi/ ReKi 
characterizes ratio of dispersion to diffusion and parame- 
ter w = —{ImAo+Trf)/ReAo, characterizes the frequency 
of the driving. 

The only difference between this equation and the 
Gin zburg-Landau equat ion for the parametric instabil- 
ity l)Coullet et all Il990(l is the coupling of the complex 
amplitude tp to the "slow mode" v which characterizes 
local dissipation in the granular layer (1/ can be inter- 
preted as coarse-grained layer's number density). This 
slow mode obeys its own dynamical equation 

dtiy ^ aV ■ {i^VliPl^) + PV^I^. (24) 

This equation describes re-distribution of the averaged 
thickness due to the diffusive flux oc — Vi', and an addi- 
tional flux cx — i/VjV'P is caused by the spatially nonuni- 
form vibrations of the granular material. Th i s cou - 
pled model was used bv lAranson and Tsimrin3 l)l998(l : 



FIG. 19 Phase diagram showing primary stable patterns de- 
rived from Eqs. . I24II . Points indicate stable oscillons 
obtained by numerical solution of Eqs. (I23II . 12411 . — a/f3, 
fj, = (v) is avera ge density, and e ~ 7 — 7e is s upercriticality 
parameter, from iTsimring and Aranson I lll997ri . 

FIG. 20 Radially-symmetric oscillon solution of 

Eqs. m, 1231 fo r 7 = 1-8. /x = 0.567 fc = = a = 

1, 77 = 5/7, from ITsimring and Aranson I il997ft . 

ITsimring and Aranson I 1)19971) to describe the pattern 
selection near the threshold of the primary bifurcation. 
The phase diagram of various patterns found in this 
model is shown in Fig. 1191 At small a{v)l3~^ (which cor- 
responds to low frequencies and thick layers), the primary 
bifurcation is subcritical and leads to the emergence of 
square patterns. For higher frequencies and/or thinner 
layers, transition is supercritical and leads to roll pat- 
terns. At intermediate frequencies stable localized solu- 
tions of Eqs. (|23|l . (12411 corresponding to isolated oscillons 
and a variety of bound states were found in agreement 
with experiment. The mechanism of oscillon stabiliza- 
tion is related to the oscillatory asymptotic behavior of 
the tails of the oscillon (see Fig. I20|l . since this under- 
lying periodic structure provides pinning for the circular 
front forming the oscillon. Without such pinning, the os- 
cillon solution could only exist at a certain unique value 
of a control parameter (e.g. 7), and would either collapse 
or expand otherwise. 

Let us note that stable localized solutions somewhat re- 
sembling oscillons have recently been found in the nonlin- 
ear Schrodinger equation wi th additional linear dissipa- 
tion and parametric driving ('Barash enkov et 

Phenomenological model (|23|l . (12411 also provides a good 
description of patterns away fro m the primary bifurca- 
tion - hexagons and interfaces ()Aranson et alV Il999a(l . 
In high-frequency limit the slow mode dynamics can be 
neglected {v becomes enslaved by ij:), and the dynam- 
ics can me described by a single parametric Ginzburg- 
Landau equation (|23|l . 

It is convenient to shift the phase of the complex order 
parameter via ip = '>pexp{i(j)) with shi2(j) — w/7. The 
equations for real and imaginary part ijj = A + iB are: 

dtA = {s-1)A-2luB-{A'^ + B^)A + V^{A-bB), 
dtB = -{s + l)B - {A^ + B'^)B + \I'^{B + bA), (25) 

where = 7^— w^. At s < 1, Eqs. (|25|l has only one triv- 
ial imiform state ^ = 0, _B = 0, At s > 1, two new uni- 
form states appear, A = ±Aq,B — Q^Aq — ^/s — 1. The 
onset of these states corresponds to the period doubling 
of the layer flight s sequ ence, observed in experiments 
()Melo et all 1 1994 [19951) and p redicted by the simple 
inela s tic ba ll model (Mehta and LuckL Il990l: ^elo et aQ, 
119941199511 . Signs ± reflect two relative phases of layer 
flights with respect to container vibrations. 
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Weakly-nonlinear analysis reveals that the uniform 
states ±Aq lose their stability with respect to finite- 
wavenumber perturbations at s < Sc, and the nonlin- 
ear interaction of growing modes leads to hexagonal pat- 
terns. The reason for this is that the non-zero base state 
A = ±Ao lacks the up-down symmetry -0 — > —-0 and the 
corresponding amplitude equations contains quadratic 
terms whi ch are known to favor hexag ons close to onset 
(see, e.g. ICross and Hohenberd l|l993|) '). In the regime 
when the uniform states A = zLAq,B = are stable, 
there is an interface solution connecting these two asymp- 
totic states. This interface may exhibit transversal in- 
stability which leads to decorated interfaces (see exper- 
imental Fig. 13?). Due to symmetry, the interfaces are 
immobile, however breaking the symmetry of driving can 
lead to interface motion. This symmetry breaking can be 
achieved by additional subharmonic driving at frequency 
//2. The interface will move depending on the relative 
phases of / and / /2 har monics of driving. This interface 
drift was predicted in ijAranson e^_a/ 1 Il999a^ and ob- 
served in the subs equent work llAranson et all Il999bfl . 
As it wa s noted bvlAranson et al\ l)l999b|) fsee also later 
work by iMoon et all ()2003(l ). moving interfaces can be 
used to separate granular material of different sizes. The 
stability and transition between flat and decorated in- 
terfaces w as studied theoretically and experimentally by 
iBlair et~aL (2000) . It was shown that non-local effects 
are responsible for the saturation of transverse insta- 
bility of interfaces. Moreover, new localized solutions 
( "superoscillons" ) were found for large accelerations. In 
contrast with conventional oscillons existing on the flat 
background oscillating with driving frequency /, i.e. in 
our notation -0 = 0, the superoscillons exist on the back- 
ground of the flat period-doubled solution 0^0. 

Another description o f the primary pattern-f orming bi- 
furcation was done by ICrawford and Riecke I (^99) in 
the framework of the generalized Swift-Hohenberg equa- 
tion 



i)V + ^^V'^ 



-eV • [(VV')^] 
(26) 



Here the (real) function ip characterizes the amplitude of 
the oscillating solution, so implicitly it is assumed that 
the whole pattern always oscillates in phase. Terms pro- 
portional to e have been added to the standard Swift- 
Hohenberg equation first introduced for description o f 
convective rolls (see, e.g. ICross and Hohenberd l)l993|) ') 
since they are known to favor square patterns, and ex- 
tended fifth-order local nonlinearity allowed to simulate 
subcritical bifurcation for R < 0. This equation also de- 
scribes both square and stripe patterns depending on the 
magnitude of e and for negative R has a stable oscillon- 
type solution. 

Even more generic appro a ch w as taken by 
IVenkataramani and Ott"1 l)l998l l200l|) who argued 
that the spatio-temporal dynamics of patterns gener- 
ated by parametric forcing can be understood in the 
framework of a discrete-time, continuous space system 



which locally exhibits a sequence of period-doubling 
bifurcations and whose spatial coupling operator selects 
a certain spatial scale. In particular they studied the 
discrete-time system 



e„+i(x) = /:[A/(6.(x))] 



(27) 



where local mapping M(^) is described by a Gaussian 
map 

M(e)=fexph(e-l)V2] 

and the linear spatial operator C has an azimuthally sym- 
metric Fourier transform 

/(fc) = sign[fc2 - exp[fc2(l - e/2kl))/2]. 

Here k is the wavenumber, kc, kg are two inverse length 
scales characterizing the spatial coupling, and f describes 
the amplitude of forcing. While this choice of the spa- 
tial operator appears rather arbitrary, it leads to a phase 
diagram on the plane {kc/ko,r) which is similar to the 
experimental one. 

Several authors JCe^^^^L 1997t 



lEggers and Rieckel Il999t iPark and MoonL l2002l) at 

tempted to develop a quasi two-dimensional fluid- 
dynamics-like continuum description of the vibrated 
sand patterns. These models deal with mass and 
momentum conservation equations which are augmented 
by specifi c constitutive relations for the mass flux and 
pressure. ICerda et al\ lfl997.'l assumed that during im- 
pact particles acquire horizontal velocities proportional 
to the gradient of local thickness, then during the flight 
that move freely with these velocities and redistribute 
mass, and during the remainder of the cycle the layer 
diffusively relaxes on the plate. The authors found that 
a flat layer is unstable with respect to square pattern 
formation, however the transition is supercritical. In 
order to account for the subcritical character of the 
primary bifurcation to square patterns, the authors pos- 
tulated the existence of a certain critical slope (related 
to the repose angle) below which the free flight initiated 
by the impact does not occur. They also observed the 
existence of localized excitations (oscillons and bound 
states), howe ver they appear e d only as transients in 
the model. iPark and MoonI l)2002j) generalized this 
model by explicitly writing the momentum conservation 
equation and introducing the equation of state for the 
hydrodynamic pressure which is proportional to the 
square of the velocity divergence. This effect provides 
saturation of the free-flight focusing instability and leads 
to a squares-to-stripes transition at highe r frequencies 
which was missing in the original model l)Cerda et al\ . 
Il997j) . By i ntroducing multiple fre e-flight times and 
contact times iPark and Moct] ||2002|) were also able to 
reproduce hexagonal patterns and superlattices. 

Full three-dimensional continuum simulations based on 
the granular hydrodynamics equations ((Sj , © , lO were 
performed by Bougie et at. (2005). Quantitative agree- 
ment was found between this description and event- 
driven molecular dynamics simulations and experiments 
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FIG. 21 Dispersion relation for stipes near the onset ac- 
cording to continuum granular hydrodynamics equations and 
molecular dynamics simulation s compared with experimental 
data, from lBougie et all i2005D . 



in terms of the wavelength dependence on the vibration 
frequency fFig l21|l ahhough the authors had to introduce 
a certain regularization procedure in the hydrodynamic 
equations in order to avoid artificial numerical instabil- 
ities for ^ — > 0. Since standard granular hydrodynam- 
ics does not take into account friction among particles, 
the simulations only yielded stripe pattern, in agreement 
with earlier molecular dynamics simulations. Further- 
more, the authors found a small but systematic difference 
(~10%) between the critical value of plate acceleration 
in fluid-dynamical and molecular dynamics simulations 
which could be attributed to the role of fluctuations near 
the onset. Proper accoimt of inter-particle friction and 
fluctuations within the full hydrodynamics description 
still remains an open problem (see more on that in Sec- 
tion ED). 

Fluctuations are expected to play a significantly 
greater role in granular hydrodynamics than in usual 
fluids, because the total number of particles involved 
in the dynamics per characteristic spatial scale of the 
problem is many orders of magnitude smaller than the 
Avogadro number. The apparatus of fluctuating hy- 
drodynamics which was developed in particular for de- 
scription o£_transi^ion__to_^olls in Ra yleigh-Benard con- 
vection ijSwift and Hohenberel 11977 ), has been r e cently 
applied to the granular patterns 1 Bougie et all l2005t 
[Goldman et al, 2004). The Swift-Hohenberg theory is 
based on the equation for the order parameter ^, 

dt^ = [e - (V2 + klf]i: -^^ + 7y(x, i), (28) 

where e is the bifurcation parameter, fco is the wavenum- 
ber corresponding to the most unstable perturbations, 
and 77 is the Gaussian 6 correlated noise term with in- 
tensity F. The Swift-Hohenberg theory predicts that 
noise offsets the bifurcation value of the control param- 
eter from the mean-field value ejv/F = to the critical 
value Cc oc F"^/^. Furthermore, the Swift-Hohenberg the- 
ory describes the transition to the linear regime which 
is expected to work far away from the bifurcation point 
for small noise intensity when the magnitude of noise- 
excited modes scales as |e — ed"^/^, while the time co- 
herence of fluctuations and the amplitude of spectral 
peaks decays as |e — ed"^. Fitting the Swift-Hohenberg 
equation (I2 8|l to match the tra n sition in vibrate d gran- 
ular laver. iBoueie et Ml (iooi); iGoldman etlJ] ll2fl04D 
found a good agreement with molecular dynamics sim- 
ulations and experiments (see, e.g., Fig l22l) . Interest- 
ingly, the magnitude of the fitted noise term in Ea. l(^ 
F « 3.5 • 10"'^ turned out to be an order of magnitude 
greater than f or convective instabili ty in a fluid near a 
critical point ijOh and Ahlersl l2003(l . This discrepancy 
could stem from the fact that the Swift-Hohenberg the- 



FIG. 22 Comparison between the Swift-Hohenberg theory 
and experiment for noise peak intensity (a) , total noise power 
(b) and the correlation time (c). Symbols - experiment, solid 
lines - Swift-Hohenberg theory, dashed lines - l inear theory 
for small noise magnitude, from IGoldman 



FIG. 23 Select patterns in a shallow fluidized bed with pe- 
riodi cally m odulated air flow for different flow parameters, 
from lLi et aL (■200a') 



ory, developed for ordinary fluids, is formally valid for 
the second-order phase transition, whereas in granular 
system the transition to square patterns is of the first 
order type. Consequently, the nonlinear terms can be 
important near the transition point and may distort the 
scaling for the noise amplitude. 

There have been attempts to connect patterns in vi- 
brated layers with the phenomenon of granular "ther- 
moconvection" . Since high-frequency vibration in many 
aspects is similar to "hot" wall, it was argued that one 
should expect granular temperature gradients, density 
inversion, and, consequently convection instability simi- 
lar to that observed in heated from below liquid layers. 
The theoretical analysis based on granular hydrodynamic 
equations lO , © , iQ supports the existence of a conve c- 
tive i n stability in a certain r a nge of parameters l)He et all . 
120021 iKhain and Meersonl |2003|) . Multiple convec- 
tion roles were observe d in molecular dynamics sim- 
ulatio ns l)Paolotti et all 12004 ISunthar and Kumaranl 
l200lj) . However, the experim ents are not conclusive 
enough l)Wildman et all l200lj) . In particular it appears 
very hard to discriminate between convection induced by 
vibration and convective flows induced by walls, se e e.g. 
l)Garcimarin et a/.L l2003 iPak and Behringerlll993|) . 

Vibrated bottom plate is not the only way t o induce 
parametric patterns in thin granular layers. iLi et all 
( 2003) demonstrated that periodically modulated airflow 
through a shallow fluidized bed also produces interesting 
patterns in the granular layer which oscillate at half the 
driving frequency (Fig. I23|l . While the physical mech- 
anism of interaction between the airflow and grains is 
quite different from the coUisional energy transfer in vi- 
brated containers, phenomenological models based on the 
principal symmetry of the problem should be able to de- 
scribe the gas-driven granular layer as well. In case of 
the parametric Ginzburg-Landau model Eq. (|23() . the 
order parameter would correspond to the amplitude of 
the subharmonic component of the surface deformation, 
and the driving term would be related to the amplitude 
of the flow modulation. Moreover, variations of the mean 
flow rate act similar to the variations of the gravitational 
acceleration in the mechanical system, which may give 
an additional means to control the state of the system. 
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FIG. 24 Stability diaRra m for avalanches in thin granular 
layers, se also Fig. |S] from lDaerr and Douadvl il999D 



VI. PATTERNS IN GRAVITY-DRIVEN DENSE 
GRANULAR FLOWS 

In this Section we overview theoretical models for var- 
ious pattern-forming instabilities in dense gravity-driven 
granular flows. 



A. Avalanches In thin granular layers 

Gravity-driven particulate flows are a common oc- 
currence in nature (dune migration, erosion/deposition 
processes, land slides, underwater gravity currents and 
coastal geomorphology) and in various industrial ap- 
plications having to do with handling granular materi- 
als, including their storage, transport, and processing. 
One of the most spectacular (and often very dangerous) 
forms of gravity-driven granular flows is the avalanche. 
Avalanches occur spontaneously when the slope of the 
granular material exceeds a certain angle (static angle of 
repose) or they can be initiated at somewhat smaller an- 
gle by applying a finite perturbation. Laboratory studies 
of avalanches are often carried out in rotating drums (see 
below) or in a chute geometry when a layer of sand is ti- 
tled at a certain fixed angle. iDaerr and Douadvl l|l999l) 
conducted experiments with a thin layer of granular mat- 
ter on sticky (velvet) inclined plane, see Fig. |^ Sur- 
prising diversity of avalanche behavior was observed in 
this seemingly simple system: triangular avalanches de- 
veloped in thin layers {h is the layer thickness) and for 
small inclination angles 0, whereas in thicker layers or 
steeper angles </> the avalanches assumed balloon shaped 
with upper edge of t he avalanche propagating up -hill, see 
Fig. El According to' RaichenbachI l)2002a '. '200?) the rear 
front of the balloon-like avalanche propagates uphill with 
the velocity roughly one half of the downhill velocity of 
the head front, and the velocity of the head is also two 
times larger than the depth-averaged flow velocity. The 
stability diagram is outlined in Fig. 1241 a granular layer 
is stable belo w solid line (so-called hstop limit according 
to jPouhaueni l)l999j) '). spontaneous avalanching was ob- 
served above the dashed line. Between dashed and solid 
lines the layer exhibits bistable behavior: finite pertur- 
bation can trigger an avalanche, otherwise the layer re- 
mains stable. The dotted line with x -symbols indicates 
the transition between triangular and balloon avalanches. 



1. Partially fluidized flows 

The avalanche dynamics described above is an exam- 
ple of a wide class of partially fluidized granular flows. In 
such flows part of grains flows past each other while other 
grains maintain static contacts with their neighbors. The 



description of such flows still represents a major challenge 
for the theory. In particular, one is faced with the prob- 
lem of constructing the constitutive relation for the stress 
tensor cr. In dense quasi-static flows a significant part of 
the stresses is transmitted through quasi-static contacts 
between particles as compared with short collisions in 
dilute flows. 

Stimulated b y the non-tr i vial avalanche dynamics i n 
experiments bvlPaerJ I200ll):lpaerr and Douadvl (|l.999), 
lAranson and Tsimrind ll200ll l2002^ suggested a generic 
continuum description of partially fluidized granular 
flows. According to this theory, the ratio of the static 
part (T** to the fluid part cr-^ of the full stress tensor is 
controlled by the order parameter p. The order param- 
eter is scaled in such a way that in granular solid p = 1 
and in well developed flow (granular liquid) p —^ 0. On 
the "microscopic level" the order parameter is defined as 
a fraction of the number of static (or persistent) contacts 
of the particles Zg to total number of the contacts Z, 
p ~ {Zs/Z) within a mesoscopic volume which is large 
with respect to the particle size but small compared with 
characteristic size of the flow. 

Due to a strong dissipation in dense granular flows the 
order parameter p is assumed to obey purely relaxational 
dynamics controlled by the Ginzburg-Landau-type equa- 
tion for the generic first order phase transition. 

Here D is the diffusion coefficient. F{p, S) is a free energy 
density which is postulated to have two local minima at 
p — 1 (solid phase) and p = (fluid phase) to account 
for the bistability near the solid-fluid transition. 

The relative stability of the two phases is controlled by 
the parameter 6 which in turn is determined by the stress 
tensor. The simplest assumption consistent with the 
Mohr-Coloumb yield criterion is to take it as a function of 
<f) = max |cr„„i/(T„„|, where the maximum is sought over 
all possible orthogonal directions m and n (we consider 
here only two-dimensional formulation of the model, an 
objective th ree dimensional g eneralization was recently 
proposed bv lGao et al\ 1)20051) 1. Furthermore, there are 
two angles which characterize the fluidization transition 
in the bulk of granular material, an internal friction an- 
gle tan^^ (j)i such that if > 0i the static equilibrium is 
unstable, and the "dynamic repose angle" tan^^ 0o such 
that at (/) < the "dynamic" phase p = 0, is unstable. 
Values of and 4>i depend on microscopic properties 
of the granular material, and in general th ey do not co- 
incide. lAranson and Tsimrin 3 l|200li |2002|) adopted the 
simple algebraic form of the control parameter (5, 

5={^^-<Pl)l{<f>\-<f^l). (30) 

Order parameter equation H29() has to be augmented by 
boundary conditions. While this is a complicated issue in 
general, a simple but meaningful choice is to take no-flux 
boundary conditions at free surfaces and smooth walls. 
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FIG. 25 Comparison of theoretical and experimental phase 
diagrams. Lines obtained from theory, sym bols depicts ex- 
perimental data from Ref. (iDaerr and Doua dv. 1999). Solid 
line and circles limit the range of existence of avalanches, 
line and triangles correspond to the linear stability bound- 
ary of the static chute, and the line and crosses denote the 
boundary between triangular and balloon avalanches. In- 
set: Schematic representation of a chute flow geometry, from 
lAranson and Tsimrind ll200ll. l2002fl 



FIG. 26 Sequence of images demonstrating the evolution 
of a triangular avalanche (a-c) and up-hill avalanche (d-f) 
obtained form numerical solution of Eqs. II32LI3H . from 
lAranson and Tsimrind (l200ll. I2002D 



model also was test ed in soft-p a rticle m o lecula r 
dynamics simulations ifVolfson et all l2003albl l2004|) . 
lOrpe and KhakhaJ ^I200,'^^ used the pa r tial fluidiza- 
tion model of f Aran son and Tsimrind. l200li l2002t 
IVolfson et aZl 12003a ') for the description of velocity pro- 
files three-dimensional shear flows in a rotating drum. 
The comparison between experimental data and theory 
shows that the partial fluidization model describes rea- 
sonably well entire velocity profile and the flow rheology, 
however experimental methods for independ ently esti- 
matin g the order parameter model are needed. iGao et all 
l)2005f) recently developed an objective (coordinate sys- 
tem independent) formulation of the partial fluidization 
theory which allows for the straightforward generaliza- 
tion to three-dimensional systems. 



and solid phase condition p — 1 near sticky or rough 
walls. 

For the flow of thin granular layers on inclined planes 
Eqs. 10, H29|l can be simplified. Using the no-slip 
boundary condition at the bottom and no-fiux condi- 
tion at the top of the layer and fixing the lowest-mode 
structure of the order parameter in the direction per- 
pendicular to the bottom of the chute (z — 0, see Inset 
to Fig. |23J), p = 1 — A{x,y)sin{'Kz/h), h is the local 
layer thickness, A(x, y) is slowly-varying function, one 
obtains equations governing the evolution of thin layer 
l| Aranson and Tsimring. 2001. 2002) : 

dth = -ad^{h^A) + ^\/ (h^AWh) , (31) 
<p 

dtA ^ XA + VlA + ^^^^A^--A'' (32) 

Sir 4 

where Wl ^ + d^, X ^ 6 ~ 1 - ir^/Ah^, a « 
0.12/z^^(7 sin (^, p, is the shear viscosity, (p is the chute 
inclination, (p = tani^. Control parameter S includes 
correction due to change in the local slope S = Sq + (3hx, 
(i « 1/(01 — 0o) ~ 1.5 — 3 depending on the value of 
(f). The last term in Eq. H31|) is also due to change of 
local slope angle and is obtained from the expansion 
if K, (p -\- hx- This term is responsible for the saturation 
of the avalanche front slope (without it the front would 
be arbitrarily stee p). While it was not inc luded in orig- 
inal publications ijAranson and Tsimriiii |2001, 20oJ, 
this term is important for large wavenumber cut-off of 
long-w ave instability observed bv lForterre and PouliauenI 
l)2003l) . see Sec. IVI.CI Numerical and analytic solutions 
of Eqs. I|31() . H32|I exhibit strong resemblance with exper- 
iment: triangular avalanches in thin layers and balloon- 
like avalanches in thicker layers, see Fig. 1201 The cor- 
responding phase diagram agrees quantitatively with an 
experimental one having only one fitting parameter (vis- 
cosity p), Fig. 123 

In subsequent work ijAranson and Tsimrinl. |2002|) . 
this theory was generalized to other dense shear 
granular flows including flows in rotating drums, 
two- and three-dimensional shear cells, etc. The 



2. Two-phase flow approach of granular avalanches 



Another approach treating near-surface granular 
flows as two-phase systems was developed by a 
numb er of authors, see e.g. llBouchaud et alV 



1994', '1995"; 'Boutrcux et a/.', ' 


1998; 


'Douadv et al\. 


Wm 


Khakhar et al. 




2001 




Meht£ 


.1994) and many 


oth- 



ers. For review on r ecent models of surface flows see 
(lAra dian et all . l2002|) . All these models distinguish 
rolling and static phases of granular flow described by 
the set of coupled equations for the evolution of thick- 
nesses of both phases, R and /t, respectiv e ly. T he phe- 
nomen o logica l theory by iBouchaud" et al\ ||1994 [1995); 
iMehtal l)l994() (often called BCRE theory) provides an 
intuitive description of the flow. In shallow granular 
layers, even simpler depth-averaged granular hydrody- 
namic equations (so-called Saint- Venant models) often 
provi d es quite accurat e descr i ption, see jPouad v et al\ 
19991 iKhakhar et all l200lt iLaieunesse et all 1200^ 
Savage and Hutted 119891 

In the most general and compact form the BCRE the- 
ory can be represented by a pair of equations for evolution 
of R and /i, 

dth = T{h,R) (33) 
dtR = VddxR-T{h,R) (34) 

where r(ft,, R) is the exchange term, or a conversion rate 
between rolling and static grains, and Vd is the downhill 
grain velocity. Physical meaning of the BCRE model is 
very simple: Eq. H33|> expresses the increase in the height 
due to deposition of roUing grains, and Eq. describes 
advection of rolling fraction by the flow with velocity 
Vd and depletion due to conversion to static fraction. 
The limitations a nd generalizati o ns of t he BCRE model 
are di scussed by lAradian et al\ l)2002() : iBoutreux et al\ 

lll998ll. 

iDouadv et al\ l)200 2f) applied the following two-phase 
model to describe avalanches in thin granular layers: 

dth + 2Udxh ^ i {tan <j) ~ p{h)) (35) 
dtC + 2Udxh = (36) 
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where C = R+h is the position of free surface, U is depth- 
averaged velocity of the flow. In addition to BCRE model 
Eqs. (|35|I . H36|I include two phenomenological functions: 
r characterizes the mean velocity gradient of a single 
bead on incline, and n{h) describes depth-dep endent fric- 
tion with the bottom. According to[Douadv et al\ l|2002ft 
a three-dimensional version of Eqs. (|35() . H36(l describes 
transition from triangular to uphill avalanches, however 
details of the transition depend sensitively on the choice 
of functions F and 

Depth- averaged d escription in th e form of Eqs. 
was used bv lBorzsonvi et al~\ l)2005() to address 
the difference between shapes of avalanches for sand and 
glass particles in a chute flow. The authors reduced Eqs. 
(|35|l . (Infill to the modifled Burgers equation 



dfh + a{h)dxh = fi{h)d^h 



(37) 



where function a{h) ^ h^^"^ and effective viscosity /i ^ 
^/h. This description connects avalanches with the 
"Burgers" shocks. Eq. (|37|l implies that all avalanches 
will eventually decay, in contrast to experiments indi- 
cating that only small avalanches decay whereas large 
avalanches grow and/or form st ationary waves (Daerr, 
l2001at iDaerr and Douadvl Il999j) . This discrepancy is 
likely due to the fact that reduction of the full model 
(|35|l . (|36|) to the single equation H37|l does not take into 
account the bistable nature of granular flows. 

While two-phase description of granular flow is sim- 
ple and rather intuitive, it can be problematic when a 
clear-cut separation between rolling and static phases is 
absent, especially near the onset of motion. The order 
parameter approach can be more appropriate in this sit- 
uation. Furthermore, the two-phase equations can be de- 
rived from the partial fluidization model described in the 
Sec. IVI.A.ll as a sharp- interface limit of the co ntinuum 
order parameter model ijAranson and Tsimrind . 12002) . 



3. Avalanche shape 

On the basis of simple kinematic considerations 
iRaichenbachI l)2002b^ suggested an analytic expression 
for the shape of triangular and balloon-like avalanches. 
For the balloon-like avalanches the shape is given by the 
envelope of the expanding circles with the center drifting 
downhill: 



y-2vt+^vT\ =[1:Vt\ , 0<t<< (38) 



where v is the velocity of the rear front. For the trian- 
gular avalanches the shape is given by the envelope of 
dilating ellipses 



vx 
2^ 



y 



3 



1 



(39) 



Here v± is perpendicular velocity. While these heuristic 
relations are consistent with experimental observation. 



FIG. 27 Top: Total area overrun by the a valanche (solid line) , 
comp ared with experimental image from ( Da err and Douadvl 



Il999f) . Bottom: superimposition of avalanche boundaries 
Riven by Eq. 11391 fo r three different moments of time, from 
IRaichenbachI j2002bD 



see Fig. 1271 their connection to continuum dynamical 
models of granular flows remains to be understood. 



B. Statistics of avalanches and sandpile model 

It is well known that in real sandpiles avalanches can 
vary widely in size. The wide distr i bution of scales in 
real avalanches stimulated iBak et al\ 1)1987 *) to introduce 
a "sandpile cellular automaton" as a paradigm model for 
the self- organized criticality, the phenomenon which oc- 
curs in slowly driven non-equilibrium spatially extended 
systems when they asymptotically reach a critical state 
characterized by a power-law distribution of event sizes. 
The set of rules which constitute the sandpile model 
is very simple. Unit size "grains" are dropped one by 
one on a one-dimensional lattice in random places and 
form vertical stacks. If a local slope (the difference be- 
tween heights of two neighboring stacks) exceeds a cer- 
tain threshold value, a grain hops from the higher to the 
lower stack. This may trigger an "avalanche" of subse- 
quent hops until the sandpile returns to the stable state. 
After that another grain is dropped and the relaxation 
process repeats. The size of an avalanche is determined 
by the number of grains set into motion by adding a single 
grain to a sandpile. This model asymptotically reaches 
a critical state in which the mean angle is equal to the 
critical slope, and avalanches have a universal power-law 
distribution of sizes, P{s) oc s~" with a ~ 1.5. 

The relevance of this model and its generalizations to 
the real avalanch es is still the ma tter of debate. The 
sandpile model bv lBak et al\ lll987|) is defined via a sin- 
gle repose angle and so its asymptotic behavior has the 
properties of the critical state for a second-order phase 
transition. Real sandpiles are characterized by two angles 
of repose and thus exhibit features of the first-order phase 
transition. Moreover, concept of self-organized criticality 
is related to a power-law distribution of avalanche sizes, 
thus reliable experimental verification of self-organized 
criticality requires accumulation of very large statistics 
of avalanche events and a large-scale experimental setup. 
Finite size effects should strongly affect the power-law 
behavior. 

Experiments with avalanches in slowl y rotating drums 
( Jaege^^aZ. '■ '198^: 'Raiche nbachl 1200 ?) and chute flows 
( Lemieux an d Durian, 2000^) do not confirm the scale- 
invariant distribution of avalanches. In other experi- 
ments with large mono-disperse glass beads dropped on a 
conical sandpile ICostello et al\ |2P0j|) claimed existence 
of the self-organized criticality with a « 1.5. Character- 
istics of the size distribution depended on the geometry 
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of the sandpile, physical and geometrical properties of 
grains, and the way the grains are dropped on the pile, 
contrary to the universal concept of self-organized criti- 
cal behavior. Self-organized criticality was also claimed 
in the avalanche statistics in three-dimensional pile of 
anisotropic grains (long rice), however a smaller scaling 
exponent a » 1.2 was measured for the avalanche size 
distribution ijAeeerter e t ai, 2004). Interestingly, rice 
piles were observed to demonstrate roughening dynam- 
ics of their surface as the distribution of active sites in 
the self-organized critical state shows a self-afHne struc- 
ture w ith the fractal exponent ds = 1.85 l)Aegerter et al\ . 
l2004l) . This is consistent with the theoretically predicted 
mapping between self-organized criticality and roughen- 
ing obser ved for example in Kardar-Parizi-Zhang model 
|Paczuski and Boettcherl [1996^) . 

One can argue that real sandpiles should not exhibit 
self-organized criticality in a strict sense due to hysteresis 
and the existence of two critical repose angles. However, 
since the difference between the angles is relatively small, 
one cannot exclude power-law type behavior in the fi- 
nite range of avalanche sizes. This circumstance possibly 
explains significant scatter in experimental results and 
scaling exponents for avalanche size distribution and the 
dependence on grain shape and material properties. 



C. Instabilities in granular chute flows 

Granular chute flows exhibit a variety of pattern- 
formin g instabilities, inc luding fingering (Mallog gi et al\ . 



2005a; Pouliauen et al 
Borffion^ 



Il997t), longitudinal vortices 
2005t Fo^^^n^Poulic^ueE , 



long surface waves IIForterre and PouliauenL 



2003^, se gregation and str atification |GravandHuttei 
'997a: Makse a/.lll997bl). etc 



iPouliauen et cH ](|l993) studied experimentally a gran- 
ular chute flow on a rough inclined plane. Experiments 
performed with polydisperse sand particles demonstrated 
fingering instability of the front propagating down the 
slope, similar to that observed in fluid films flowing down 
inclined plane ijTroian et alV ^989; Zhou et al., 2005). 
However, similar experiments with smooth monodisperse 
glass beads exhibited no instability. The authors argued 
that the instability was due to a flow-induced size segre- 
gation in a polydisperse granular matter. The segrega- 
tion indeed was foun d near the avalanche front . However, 
similar experiments ijMalloggi et alV l2005albr) showed a 
fingering front instability without a significant size segre- 
gation. Thus, the question of the mechanism of fingering 
instability is still open. 

Experimen ts by |B6rzs6nvi and Ecke| l)2005j) ; 
iRirterre and PouliauenI HoOl, .2002) show the de- 
velopment of longitudinal vortices in rapid chute flows, 
see Fig. |H| The vortices develop for large inclination an- 
gles and large flow rates in the regime of accelerating flow 
when the flow thickness decrease s and the mean flow ve- 
locity increases along the chute. iForterre and PouliauenI 



FIG. 28 Density profiles 1^(2) as function of distance form the 
chute bottom z fo r different values of mean flow velocity, from 
IForterre and Pouhauen. C200^ 



FIG. 29 Phase diagram in mean density [v) and flow thick- 
ness (ft) plane delineating different fl ow instabilities. Smaller 
nu co rresponds to faster flow, from IForterre and PouliauenI 
(120021) 



l|2001|) proposed an explanation of this phenomenon in 
terms of granular "thermoconvection" . Namely, rapid 
granular flow has a high shear near the rough bottom 
which leads to the local increase of granular temperature 
and consequently creates a density inversion. In turn, 
the density inversion trigger a convection instability 
similar to that in ordinary fluids. The critical instability 
wavelength Ac is determined by the depth of the layer h 
(in experiment Ac ~ 3/i). 

In a subsequent work IForterre and PouliaueiJ l)2002j) 
studied the formation of longitudinal vortices and the 
stability of granular chute flows in the framework of gran- 
ular hydrodynamics Eqs.lQ-Q. The inverse density pro- 
file appears when a heuristic boundary condition at the 
bottom relating slip velocity and heat flux is introduced. 
Steady-state solution of Eqs.|(Sl)-(IZ|) indeed yields an in- 
verse density profile (Fig. I28|) which turns out to be 
unstable with respect to short-wavelength perturbations 
for large flow velocities, see Fig. j^Hl While the linear sta- 
bility analysis captured many important features of the 
phenomenon, there are still open questions. The stabil- 
ity analysis was performed for the steady flow whereas 
the instability occurs in the regime of accelerating flow. 
Possibly due to this assumption the linear stability anal- 
ysis yielded oscillatory instability near the onset of vor- 
tices, whereas for the most part, the vortices appear to 
be steady. Another factor which is ignored in the theory 
is the air drag. The high flow velocity in the experiment 
(about 1-2 m/sec) is of the order of the terminal velocity 
of an individual grain in air, and therefore air drag may 
affect the granular flow. 

IForterre and PouliauenI (|2003) presented an experi- 
mental study of the long-surface-wave instability devel- 
oping in granular flows on a rough inclined plane. Fig. 
1501 This instabi li ty was known from previous stud- 
ies l|Davie^. '19901: ISavagel ^7^, however no precise 
characte rizatio n of the i nstab ility had been performed. 
IForterre and PouliauenI l)2003|) measured the threshold 
and the dispersion relation of the instability by impos- 
ing a controlled perturbation at the entrance of the flow 
and measuring its evolution down the slope, see Fig. 
The results are compared with the prediction of a linear 



FIG. 30 Long-surface wave instability observed in flow of sand 
down rough incline, from iFortcrrc and Pouliauen. L2003.) 
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FIG. 31 Experimental dispersion relation for the long surface 
wave instability. Shown spatial gro wth rate as function of 
the fr equency of forcing wave, from iForterre and PouliouenI 



FIG. 32 The growth rate of small perturbations a vs wave- 
length k derived from Eqs. for /3 = 2, a = 
0.025,5 = 1.1. Instabihty occurs near hsto-p curve in Figs. 
I24I25I (h — 2.9) and disappears with further increase of h. 

stability analysis conducted in the framework of depth- 
averaged Saint- Venant-type equations similar to those 
described in Sec. IVI.A.2I 

dth + d,{uh) = (40) 
dt{uh) + adxiu^h) = (tan0 — fi{u,h) — dxh ) ghcos{9) 

where h is local thickness, 9 is inclination angle, u 
is depth-averaged flow velocity, ^{h, u) is a function 
describing effective depth and velocity dependent bot- 
tom friction, a ~ 0{1) is a constant determined by 
the velocity profile wi t hin t he layer. According to 
IForterre and P ouliauenI l)200,'^ . the instability is similar 
to the long-wave instability observed in classical fluids 
but with characteristics that can dramatically differ due 
to the specificity of the granular rheology. The theory is 
able to predict quantitatively the stability threshold and 
the phase velocity of the waves but fails to describe the 
observed cutoff of the instability at high wavenumbers. 
Most likely, one needs to include higher order terms, such 
as d^h in the first Eq. H40|l in order to account for the 
cutoff. 

The order parameter theory based on Eqs. (|31|l . (|32() 
also reproduces the long-surface wave instability. Fur- 
thermore, linearizing Eqs. (|31|l . (|32|) near the steady fiow- 
ing solution A = Aq + cl exp [at + ikx],h ~ Hq + H exp [at + 
ikx], after simple algebra one obtains that the growth 
rate of linear perturbations a is positive only in a band 
restricted by some critical wavenumber and only in the 
vicinity of hstop, see Fig. (221 With the increase of h, i.e. 
the granular flux, the instability disappears, in agree- 
ment with experiments. The nonlinear saturation of the 
instability results in the development of a sequence of 
avalanches, which is generally non-periodic, see Fig. 
The structure shows slow coarsening due to merging of 
the avalanches. This instability is a possible candidate 
mechanism of the formation of inhomogeneous deposit 
structure behind the front of an avalanche. 

Conwav et al. (2003) studied free-surface waves in 
granular chute flows near a frictional boundary. The 



FIG. 33 Typical profiles of hight h and order parameter A 
in the regime of long-surface wave instability for (5 = 2, 
a = 0.025, 5 = 1.1. Starting from generic initial conditions 
h = ho, A = const plus small noise, a sequence of avalanches 
develops. 



experiments showed that the sub-boundary circulation 
driven by the velocity gradient plays an important role 
in the pattern formation, suggesting a similarity between 
wave generation in granular and fluid flows. 

A Kelvin-Helmholt z-like shear ins t ability in chute 
flows was observed bv lGoldfarbs et o/l l)2002() . when two 
streams of sand flowing on an inclined plane with dif- 
ferent velocities were in side-by-side contact with each 
other. For sufficiently high chute angles and shear 
rates the interface remains flat. The instability of the 
interface develops when the chute angle and/or the 
shear rate is reduced. This instability has been repro- 
duced in soft-particle molecular dynamics simulations by 
ICiamarra et al\ l|2005|) who also observed that in a poly- 
disperse medium this instability leads to grain segrega- 
tion (See below Sect. rvll|) . 



D. Pattern-forming instabilities in rotating cylinders 

Granular media in rotating horizontal cylinders 
(drums) often show behavior similar to chute flows. For 
very small rotation rates (as defined by small Froude 
number Fr = uP'R/ g, where uj is the angular velocity of 
drum rotation and R its radius), well separated in time 
avalanches occur when the slope of the free surface ex- 
ceeds a certain critical angle 9c whereby di minishing this 
angle to a smaller static repose angle 9c llJaeeer et all 
'l98<j^ 'R,aichenbachLll990HTegzes 6^ oil 12002. '2003V The 
difference between 9c and 9s is usually a few degrees. At 
an intermediate rotation speed, a continuous fiow of sand 
emerges instead of discrete avalanches through a hys- 
teretic transition, similar to the tr ansition in chute flows 
at lar ge rates of grain deposition ijLemieux and DurianL 
120001) . In the bulk, the granular material rotates almost 
as a solid body with some internal slipping. As mov- 
ing grains reach the free surface they slid e down within 
a thin near-surface layer ijZik et all (l994|) (see sketch in 
Fig. |5J. The surface has a nearly flat shape; the arct- 
angent of its average slope deflnes the so-called dynamic 
angle of repose 9d- 

There are various models addressing the nature of 
the transition from discrete avalanches to the continuum 
flow. L inz and Han ggi ( 1995) proposed a phenomenolog- 
ical model based on a system of equation for the angle of 
repose (j) and mean flow velocity v 

V — g (sin0 — fc(w) COS0)) x(0, ■f^) 

(j) = UJ — av (41) 

where Hi is the rotation frequency of the drum, k{v) — 
bo + h2V^ is the velocity dependent friction coefficient, 
x((/), v) is some cut-off function, and a, bo, 62 are parame- 
ters of the model. Despite the simplicity, the model yields 
qualitatively correct transition from discrete avalanches 
to continuous flow with the increase of rotation rate uj, 
and also predicts logarithmic relaxation of the free sur- 
face angle in the presence of vibration. 
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The transition from avalanches to flow naturally arises 
in the framework of t he p artial fluidization theory, 
l|Aranson and Tsimrinei |2nn2) • In this case one can de- 
rive a system of coupled equations for the parameter S 
(which is related to the surface local angle 0, see Eq. 
()30|l ) and the width of fluidized layer zq, 

dtzo = d^zo + F{zo, S) - vdsZo 
dtS = uj + dlJ (42) 

where s is the coordinate along the slope of the granular 
surface inside the drum, J = /(zq) is the downhill flux of 
grains, v is averaged velocity in flowing layer, and func- 
tions F, / and vq are derived from Eq. (|29|l . This model 
bears resemblance to the BCRE-type models of surface 
granular flows which were app lied to rotating drums by 
iKhakhar el ^ lll99/D : lMaks?lll999D . 

Eqs. exhibit stick-slip type oscillations of the sur- 
face angle for slow rotation rates and a hysteretic transi- 
tion to a steady flow for larger rates. Eqs. (|42ll yields the 
following scaling for the width of the flowing layer in 
the middle of the drum vs rotation fr equency: zp 
which is consistent with experiment l*Te gzes et alV l2002l 
After integration over s Eqs. (|42|l can be reduced 
to a system of two coupled equations for averaged drum 
angle {S) and avera ged flow thickness ( zp) so mewhat sim- 
ilar to the model of|Linz and Hanggl l)l995l) . 

Granular flows in long rotating drums u nder certain 
conditions a lso ex hibit fingering instability ijFried et alV . 
[1998; Sh(3, 12002^ . Similarity between fingering m ro- 
tating drums and chute fiows ijForterre and Pouliauenl 
1^002) suggests that mechanisms described in the Section 
IVl.CI can be responsible for this effect, see also Section 

IVITl 

VII. MODELS OF GRANULAR SEGREGATION 

One of the most fascinating features of heterogeneous 
(i.e., consisting of different distinct components) granu- 
lar materials is their tendency to segregate under external 
agitation rather than to mix, as one would expect from 
the naive entropy considera tion. This pro perty is ubiq- 
uitous in Nature (see, e.g. l)lversonl ll997l)l and has im - 
portant technological implications ijCooke et alV Il976j) . 
In fact, some aspects of segregation of small and large 
particles can be understood on equilibrium thermody- 
namics grounds l)Asakura and OosawaL Il958() . Since the 
excluded volume for small particles around large ones 
becomes smaller when large grains clump together, sep- 
arated state possesses lower entropy. However, granu- 
lar systems are driven and strongly dissipative, and this 
simple equilibrium argument can only be applied qual- 
itatively. The granular segregation is more widespread 
than it would be dictated by thermodynamics. In fact, 
any variation in mechanical properties of particles (size, 
shape, density, surface roughness, etc.) may lead to their 
segregation. At least for bi-disperse rapid dilute fiows 
the granular segregation can be rigorously treated in the 



FIG. 34 Granular s tratification in a flow down heap, from 
iMakse et ai\ Jl997bD . 

framework of kinet ic theory o f dissipative gases, see Sub- 
sec. IIILAI iJenkins and Yoon ( 2Q02i) employed kinetic 
theory for a binary mixture for spheres or disks in grav- 
ity and derived a simple segregation criterion based on 
the difference of partial pressures for each type of parti- 
cles due to the difference in size and/or mass. 

Segregation has been observed in most flows of 

S yanular mixtures, including gra nular convection 
Knight alV Il993t) . hopper flows ( Gr ay and Hutteit 
1997at iMakse et alV Il997bt JSamadani et all Il999t 
Samadani and KudrolliL l200l|), flo ws in rotating 
drums (IChoo et alV Il997t iHill 119971: IZik et alV [1994^, 
and even in freely cooli ng binary granular gases 
ijCatuto and Marconi l2004|l . Segregation among large 
and small particles due to shaking has been termed 
"Brazfl nut effect" l|Rosato et ad Il987|) . The phe- 
nomenon of granular segregation was discovered long 
time ago, and several "microscopic" mechanisms have 
been proposed to explain its nature, inclu ding inter- 
parti cle coUisions l|Brownl fl939fl . percolation l)Williamsl 
Il976j) . and others. In certain cases, separation of grains 
produces interesting patterns. For example, if a binary 
mixture of particles which differ both in size and in 
shape is poured down on a plane, a heap which consists 
of thin alternating layers of separated particles is formed 
l)Grav and Hutteil. Il997at iMakse ei"an . Il997bl) . see Fig. 
1341 Rotating of mixtures of grains with different sizes in 
long drums produce s well separated b ands o f pure mono- 
dispe r se p articles dChicarro et al\ Il997l: IChoo et'ah . 
119971: IffiiiL 119971: IZik et all 11991 . Fig. [TTl In this 
Section we only address models of pattern formation 
due to segregation (stratification and banding), without 
discussing other manifestations of granular segregation. 



A. Granular stratification 

Granular stratification occurs when a binary mix- 
ture of particles with differ ent physical prop e rties is 
slowly poured on a plate l|G ray_aiid Hutted . Il997at 
iKoeppe et all . [1998: Ma kse et a^J . ll997b() . More specifi- 
cally, it occurs when larger grains have additionally larger 
roughness resulting in a larger repose angle and the fiux 
of falling particles is small enough to cause intermittent 
avalanches down the slopes of the heap. The basic mecha- 
nism of strati fication is related to t he avalanches acting as 
kinetic sieves l)SavagellT988lll993(l . During an avalanche, 
voids are continuously being created within fiowing near- 
surface layer, and small particles are more likely to fall 
into them. This creates a downward fiux of smaller par- 
ticles which is compensated by the upward fiux of larger 
particles in order to maintain a zero total particle fiux 
across the flowing layer. Other models of granular seg- 
regation in a thin flowing layer ijPolgunin et all Il998l: 
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FIG. 35 Cellular automa ta model of granular stratification, 
from IMakse et al\ Jl997aD . 



iKhakhar "^119971 ITool lead to a similar result. Each 
avalanche leads to the formation of a new pair of layers 
in which the grains of different sorts are separated (see 
Fig. 1341) . This pair of layers grows from the bottom of 
the pile by upward propagation of a kink at which small 
particles are stopped underneath large ones. However, 
when the larger particles were smooth and small parti- 
cles were rough, instead of stratification only large scale 
segregation with small particles near the top and large 
p articles near the bottom was observed. 

iMakse et al\ l|l997allJ ) proposed a cellular automata 
model whi c h gen eralized the classical sandpile model 
l)Bak et all , Il987^ (see Section rVLB|) . In this model, a 
sandpile is built on a lattice, and rectangular grain have 
identical horizontal size but different heights (see Fig. 
135^1. Grains are released at the top of the heap sequen- 
tially, and they are allowed to roll down the slope. A 
particle would become rolling if the local slope (defined 
as the height difference between neighboring columns) 
exceeds the repose angle. To account for difference in 
grain properties, four different repose angles 0af3 were in- 
troduced for grains of type a rolling on a substrate of 
type f3 {a,P G {1;2} where 1 and 2 stand for small and 
large grains, respectively). Normally, 621 < 6*12 because 
of the geometry (small grains tend to get trapped by 
large grains), and one-component repose angles usually 
lie within this range, 621 < On, 622 < 612- However the 
ratio of 9ii, 622 depends on the relative roughness of the 
grains. For 621 < 9ii < 622 < 6*12 (large grains are more 
rough), the model yields stratification in agreement with 
experiment (Fig. I^Sb). If, on the other hand, ^22 < ^11 
(which corresponds to smaller grains being more rough) , 
the model yields only large-scale segregation: large par- 
ticles collect at the bottom of the sandpile. 

This physical model can also be recast in the form 
of continuum equations (B outreux and de Genneslll996t 
iMakse et 

iO, ri997 a) which generalize the single-species 
BCR E model of surface granular flows l)Bouchaud et all 
1X9941) (see Section Ivl: 



dtRa = 

dth = 



-VadxRa + Fo 



(43) 
(44) 



where Ra{x, t), Va are the thickness and velocity of rolling 
grains of type a, h{x, t) is the instantaneous profile of 
the sandpile, and Fq, characterizes interaction between 
the rolling grains and the substrate of static grains. In 
the same spirit as in the discrete model, the interaction 
function F^ is chosen in the form 
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(45) 



FIG. 36 Granular avalanche-induced stratifi cation in rotating 
drum observed for low rotation rates, from lGrav and Hutted 
tl997ay 



form of the interaction terms implies that the grains of 
type a become rolling if the local slope exceeds the repose 
angle 6a{4>p) for this type on a surface with composition 
(f)l3{x,t). Assuming that the generalized repose angles 
9a{4>p) are linear functions of the concentration 
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(46) 
(47) 



Eqs. (|43|I - H45|) possess a stationary solution in which 
the heap is separated into two regions where ^2(02) < 
9 < 6*1 (^2) and 9 < 6*2(^2) < 6*1 (02)- This solution 
corresponds to small grains localized near the top and 
small grains near the bottom with a cont inuous tran- 
sition b etween the two regions. However, Makse et all 
l)l997bj) showed that this stationary solution is unstable 
if (5 = 022 — &11 > and gives rise to the stratification 
pattern. 

Similar effect of stratification patterns was observed 
experimentally in a thin slowly rotating drum which 
is m ore than half fill e d with a similar binary mix- 
ture l|Grav and Hutteil Il997a(l . see Fig. EHI Periodic 
avalanches, occurring in the drum, lead to formation of 
strata by the same mechanism described above. 



B. Axial segregation in rotating drums 

The most common system in which granular segrega- 
tion is studied is a rotating drum, or a partially filled 
cylinder rotating around its horizontal axis (see Section 
IVI.Pp . When a polydisperse mixture of grains is ro- 
tated in a drum, strong radial segregation usually occurs 
within just a few revolutions. Small and rough parti- 
cles aggregate to the center {core) of the drum, large and 
smooth particles rotate around the core (see Figs. [TUl 
and|51). Since there is almost no shear flow in the bulk, 
the segregation predominantly occurs within a thin flu- 
idized near-surface layer. For long narrow drums with 
the length much exceeding the radius, radial segregation 
is often followed by the axial segregation occurring at 
later stages (after several hundred revolutions) when the 
angle of repose of small particles exceeds that of large 
particles. As a result of axial segregati on, a pattern of 
well segregated bands is formed (^HilL il997t IZik et all . 



Here (j)a{x,t) is the volume fraction of grains of type a, 
and 9i = —dxh is the local slope of the sandpile. This 



Il994j) (see, e.g.. Fig. [TT|) which slowly merge and coarsen. 
Depending on the rotation speed, coarsening can either 
saturate at a certain finite bandwidth at low rotation 
speed s when discrete aval a nches provide granular trans- 
port l|Frette and Stavan^ Il997|) or at higher rotation 
rates in a continuous flow regime it can lead to a fi- 
nal state in which all sand is separated in two bands 
l|Fiodor and Ottind.l2005lZik a/Ill994j) . 
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The axial segregation has been weh known in the en- 
gineering communit y, it was appa rently first observed 
by Oyama in 1939 ijOvamaL Il939l) . The mechanism of 
axial segregation is apparently related to the different 
friction properties of grains which lead to different dy- 
namical angles of repose. The latter are defined as the 
angle of the slope in the drum corresponding to contin- 
uous flow regime, however in real drums the free surface 
has a more complicated S-sh apc (Elpcrin and Vikhansky, 
1998t iMaksel Il999t lOrpe and Khakhaa .200L .Zik et aU 
1994 ). According to IZik et al\ lll994|) (see also l|Levin'3. 
1999|)), if there is a local increase in concentration of par- 
ticles with higher dynamic repose angle, the local slope 
there will be higher, and that will lead to a local bump 
near the top of the free surface and a dip near the bot- 
tom. As the particles tend to slide along the steepest 
descent path, more particles with higher repose angle 
will accumulate in this location, and the instability will 
develop. Zik et al. (1994) proposed a quantitative con- 
tinuum model of axial segregation based on the equa- 
tion for the conservation equation for the relative con- 
centration of the two components ("glass" and "sand"), 
c(z,t) = {pA - Pb)/{pa + Pb), 

dtC = - — {tan0A-t&i-iOB)d,il-c^){{l + yl)^). (48) 

PT Vx 

Here x and z are Cartesian horizontal coordinates across 
and along the axis of the drum, y{x, z,t) describes the in- 
stantaneous free surface inside the drum, px = Pa + Pb, 
C is a constant related to gravity and effective viscosity 
of granular material in the flowing layer. The term in 
angular brackets denotes the axial flux of the glass beads 
averaged over the cross-section of the drum. The pro- 
file of the free surface in turn should depend on c{z,t). 
If ((1 + yDyc/Vx) < 0, hnearization of Eq.(gHJl leads to 
the diffusion equation with negative diffusion coefficient 
which exhibits segregation instability with growth rate 
proportional to the square of the wavenumber. It is easy 
to see that the term in angular brackets vanishes for a 
straight profile yx — const{x). However, for the exper- 
imentally observed S-shaped profile of the free surface 
Izik et al. ( 1994J calculated that the instability condition 
is satisfied when the drum is more than half full. While 
experiments show that axial segregation in fact observed 
even for less than 50% filling ratio, the model gives a good 
intuitive picture for the me chanism of t he instab ili ty. 

Recent experiments (IChoo et ali Il997l Il998t 
iFiodor and Ottind. 120031: iHill L Il997|) have re vealed 
interesting new features of axial segregation. 'Hill' (1997*) 
perfo rmed magnetic resonance imaging studies (Hill. 
Il997^ which demonstrated that in fact the bands of 
larger particles usually have a core of smaller particles. 
More recent experiments by iFjodor and Ottino (2003) 
showed that small particles formed a shish kebab-like 
structure with bands connected by a rod-lik e core, while 
large partic les formed disconnected rings. IChoo et al\ 
1119971 f 1998(1 found that at early stages, the small-scale 
perturbations propagate across the drum in both direc- 



FIG. 37 Space-time diagram of the surface of long rotating 
drum demonstrating oscillatory size segregation. The plot 
shows the full length of the drum and extends over 2,400 sec, 
or 1,850 revolutions. Black bands correspond to 45-250 fj,m 
black san d and white bands correspond to 300-850 /im table 
salt, from lChoo et al\ (|l993). 

tions (this was clearly evidenced by the ex periments on 
the dynamics of pre-segregated mixtures (|Choo et al\ . 
119971) '!. while at later times more long-scale static 
perturbations take over and lead to the emergence of 
quasi-stationary bands of separated grains (see Fig. 1^ . 
The slow coarsening pr ocess can be a ccelerated in a 
drum of a helical shape l|Zik et alV Il994|) . Alternatively, 
the bands can be locked in an axisym metrical drum w ith 
the radius modulated along the axis l|Zik et adll99^ . 

In order to account for the oscillat ory behavior of ax- 
ial seg r egation at the in i tial sta ge , Aranson and Tsimrind 
(ll999ii:lAransoii et al\ (|l999bj) generalized the model of 
IZik et al\ lll994|) . The key assumption was that be- 
sides the concentration difference, there is an addi- 
tional slow variable wh i ch is in volved in the d y namics . 
lAranson and Tsimrind l)l999() : lAranson et al\ (|l999b|) 
conjectured that this variable is the instantaneous slope 
of the granular material (dynamic angle of repose) which 
unlike Ea. H48|) is not slaved to the relative concentration 
c, but obeys its own dynamics. The equations of the 
model read 

dtc = -d,{-Dd,c + g{c)d,9), (49) 

dtO = a{n-9 + f{c)) + Dgd,,9 + -rd,,c. (50) 

The first term in the r.h.s. of Ea. H49|l describes diffusion 
flux (mixing), and the second term describes differential 
flux of particles due to the gradient of the dynamic repose 
angle. This term is equivalent to the r.h.s. of Ea. H48|) 
with a particular function g(c) = Go{l ~ (?)■ For sim- 
plicity, the constant Gq can be eliminated by rescaling of 
distance x x/y/Go- The sign + before this term means 
that the particles with the larger static repose angle are 
driven towards greater dynamic repose angle. This differ- 
ential flux gives rise to the segregation instability. Since 
this segregation flux vanishes with g{c) \c\ — > 1 (which 
correspond to pure A or B states), it provides a natural 
saturation mechanism for the segregation instability. 

Parameter in the second equation is the normalized 
angular velocity of the drum rotation, and /(c) is the 
static angle of repose whi ch is an increa sing function of 
the relative concentration l|Koeppe et a/.l. ll9981 (for sim- 
plicity it can be assumed linear, /(c) = F+foc). The con- 
stant F can be eliminated by the substitution 6 — F . 
The flrst term in the r.h.s. of Ea. H50|l describes the lo- 
cal dynamics of the repose angle (fi increases the angle, 
and —0 -\- f{c) describes the equilibrating effect of the 
surface ffow), and the term DgdxxO describes axial dif- 
fusive relaxation. The last term, jdxxC, represents the 
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FIG. 38 Dispersion relation \{k) for segregation instability 
(left) and comparison of the frequen cy of band oscillations 
Im(X) with experiment (right), from lAranson and Tsimrind 
1119981) . 



FIG. 39 Space-time diagrams demonstrating initial 
band oscillations and c onsequent coarsening, from 
lAranson and Tsimrind (119991) . 



lowest-order non-local contribution from an inhonioge- 
neous distribution of c (the first derivative dxC cannot 
be present due to reflection symmetry x — > —x). This 
term gives rise to the transient oscillatory dynamics of 
the binary mixture. 

Linear stability analysis of a homogeneous state c = 
coi^o = + /qCo reveals that for go/o > aZ? long-wave 
perturbations are unstable, and if goj > {Dq — 13)^/4, 
short-wave perturbations oscillate and decay (two eigen- 
values Ai^2 are complex conjugate with negative real 
part), see Fig. |2HI Thi s agrees with the general phe- 
nomenology observed bv lChoo et al\ ()l997(l both quali- 
tatively and even quantitatively (Fig. 138b). The results 
of direct numerical solution of the full model (|49|l . (|50() 
are illustrated by Fig. 123 It shows that short-wave 
initial perturbations decay and give rise to more long- 
wave non-oscillatory modulation of concentration which 
eventually leads to well-separated bands. At long times 
(Fig. |2nt>) bands exhibit slow coarsening with the 
number of bands decreasing logar it hmically with time 
(see a ls o iFiodor and OttiiioT l)2003(l ; iFrette and Stavan^ 
l|l997|) : iLevitanI |)1998() ). This scaling follows from the 
exponentially weak inter action between interfac e s sep- 
arating different ban ds ijAranson and Tsimrind . Il999t 
iFraerman et aLUl997j) . 

While these continuum models of axial segregation 
showed a good qualitative agreement with the data, re- 
cent experimental observations demonstrate that the the- 
oretical u nderstanding of axial segre gation is far from 
complete ijOttino and K hakharl l2000(l . The interpreta- 
tion of the second slow variable as the local dynamic 
angle of repose implies that in the unstable mode the 
slope and concentration modulation should be in phase, 
whereas in the decaying oscillatory mode, these two 
fields have to be s hifted in phase. Further experiments 
l)Khan et al\ . l2004j) showed that while the in-phase rela- 
tionship in the asymptotic regime holds true, the quadra- 
ture phase shift in the transient oscillatory regime is not 
observed. That lead Khan et ai, (2004) to hypothesize 
that some other slow variable other than the angle of 
repose (possibly related to the core dynamics) may be 
involved in the transient dynamics. However, so far ex- 
periments failed to identify which second dynamical field 
is necessary for oscillatory transient dynamics, so it re- 
mains an open problem. Another rece nt experimental 
observation by iKhan and Morri3 l|2005l) suggested that 
instead of normal diffusion assumed in Eas. (|49|l . (|50|l . a 



slower subdiffusion of particles in the core takes place, 
(r) ~ t"' with the scaling exponent 7 close to 0.3. The 
most plausible explanation is that the apparent subd- 
iffusive behavior is in fact a manifestation of nonlinear 
concentration diffusion which can be described by equa- 
tion 

dtc = d,D{c)d,c. (51) 

For example, for the generic concentration-dependent dif- 
fusion coefficient D ^ the asymptotic scaling behavior 
of the concentration c{z,t) is given by the self-similar 
function c ^ F{z/t") /t" for t ^ 00 with the scaling ex- 
ponent a = 1/3 close to 0.3 observed experimentally. Ex- 
perimentally observed scaling function F{x/t°') appears 
to be consistent with that of Eq. H51I) except for the 
tails of the distribution where c ^ and the assumption 
D ^ c is possibly violated. Normal diffusion behavior 
corresponding to D = const and a = 1/2 is in strong 
disagreement with the experiment. 

iNewev et all l)2004j) conducted studies of axial segre- 
gation in ternary mixtures of granular materials. It was 
found that for certain conditions bands of ternary mix- 
tures oscilla te axially. In contrast to the experiments of 
IChoo et all (1997, 1998), the oscillations of bands appear 
spontaneously from initially mixed state, which strongly 
indicates the supercritical oscillatory instability. While 
in binary mixtures the oscillations have the form of peri- 
odic mixing/demixing of bands, in the ternary mixtures 
the oscillations are in the form of periodic band displace- 
ments. It is likely that the mechanism of band oscillations 
in ternary mixtures is very different from that of binary 
mixtures. One of possible explanations could be that the 
third mixture component provides an additional degree of 
freedom necessary for oscillations. To demonstrate that 
we write phenomenological equations for the concentra- 
tion differences Ca = ci — C2 and Cb — C2 — c^, where 
ci,2,3 are the individual concentrations. By analogy with 
Eq. we write the system of coupled equations for the 
concentration differences Ca.b linearized near the fully 
mixed state: 

dtCA = DAdlCA+liAdlCB. 

dtCB = DbO'^^Cb + f^Bd^CA- (52) 

If the cross-diffusion terms have opposite signs, i.e. 
fiAfJ'B < 0, the concentrations Ca.b will exhibit oscil- 
lations in time and in space. Obviously this mechanism 
is intrinsic to ternary systems and has no counterpart in 
binary mixtures. 

Parallel to the theoretical studies, molecular dynam- 
ics simulations have been performed fRapap orl l2002t 
IShoichi, 1998; Taberlct et ai, 2004). Simulations al- 
lowed researchers to probe the role of material param- 
eters which would be diffic ult to acc e ss in laboratory 
experiments. In particular, lRapapor'3 (|2002) addressed 
the role of particle-particle and wall-particle friction co- 
efficients separately. It was found that the main role 
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FIG. 40 Granular Taylor vortices observed in vertically ro- 
tating air-fluidized cylinder filled with binary mixture. Left 
image depicts entire cylinder height and width, and right im- 
age shows the dependence of concentration of small particles 
along the bed height, from IShinbrod l(2C 



is played by the friction coefScients between the parti- 
cles and the cylinder walls: if the friction coefficient be- 
tween large particles and the wall is greater than that 
for smaller particles, the axial segregation always oc- 
cur irrespective of the ratio of particle-particle friction 
coefficients. However, if the particle-wall coefficients 
are equal, the segregation may still occur if the friction 
amo ng large particles is gr eater than among small parti- 
cles. iTaberlet et oil l|2004j) studied axial segregation in a 
system of grains made of identical material differing only 
by size. The simulations revealed rapid oscillatory mo- 
tion of bands, which is not necessarily related to the slow 
b and appearence/disappea r ence observed in exp eriments 
of lChoo fi^ o.l\ lll997t 119981 1'Fio dor and Ottinol ll200.'^ . 

A different type of discrete el ement modelling of ax- 
ial segregation was proposed by lYanagitJ l)l999l) . This 
model builds upon the lattice-based sandpile model and 
replaces a rotating drum by a three-dimensional square 
lattice. Drum rotation is modelled by correlated displace- 
ment of particles on the lattice: particles in the back are 
shifted upward by one position, and the particles at the 
bottom are shifted to fill the voids. This displacement 
steepens the slope of the free surface, and once it reaches 
a critical value, particles slide dowi i according to the rules 
similar to the sandpile model of iBak et all |)1987|) but 
taking into account different critical slopes for different 
particles. This model despite its simplicity reproduced 
both radial and axial segregation patterns and therefore 
elucidated the critical components needed for adequate 
description of the phenomenon. 



C. Other examples of granular segregation 

As we have seen in the previous Section, granular seg- 
regation occurs in near-surface shear granular flows, such 
as in silos, hoppers, and rotating drums. However, other 
types of shear granular flows may also lead to segrega- 
tion. For example, Taylor-Couette flow of granular mix- 
tures between two rotating cylinders leads to formation of 
Taylor vortices and then in turn to segregation patterns 
|ShinbrotV2004'), sec Fig. EOI 

.Pouhguen et al. (199.7.) observed granular segregation 
in a thin granular flow on an inclined plane. In this case, 
segregation apparently occurs as a result of an instability 
in which concentration mode is coupled with hydrody- 
namic mode. As a result, segregation occurs simultane- 
ously with a fingering instability of the avalanche front 
(Fig. [7|l . As an implicit evidence of this relation between 
segreg ation and fingering instability, iPouliauen et al\ 
lll997li found that mono-disperse granular material does 



not exhibi t fingering instability. However, other exper- 
iments ( S henL |2002|) indicate that in other conditions 
(more rapid flows), flngering instability may occur even 
in flows of mono-disperse granular materials. Thus, the 
segregation is likely a consequence rather than the pri- 
mary cause of the fingering instability. 

An interesting recent example of pattern formation 
caused by granular segregation in a horizontally shaken 
layer o f binar y granular mixtu re was presented bv lMuUinl 
l|2000ll2002|) : lReis and MullinI l)2002|) . After several min- 
utes of horizontal shaking with frequency 12.5 Hz and 
displacement amplitude 1 mm (which corresponds to the 
acceleration amplitude normalized by gravity F = 0.66), 
stripes were formed orthogonal to the direction of shak- 
ing. The width of the stripes was growing continuously 
with time as o? oc i"'^^, thus indicating slow coarsening 
(Fig. |SJ . This power law is consistent w ith the diffusion- 
mediated mechanism of stripe merging. iR.eis and M^iUinl 
(,2002) argued on the basis of experimental results on pat- 
terned segregation in horizontally shaken layers that the 
segregation bears features of the second-order phase tran- 
sition. Critical slow-down was observed near the onset of 
segregation. The order parameter is associated with the 
combined filling fraction C, or the layer compacity. 



C = 



S 



(53) 



where Ns^i are numbers of particles in each species, As^i 
are projected two-dimensional areas of th e respective in- 
dividu al particles, and S is the tray area. lEhrhard et al\ 
l|2005l) proposed a simple numerical model to describe 
this phenomenon of segregation in horizontally vibrated 
layers. The model is based on a two-dimensional system 
of hard disks of mass and radius Ra (a = 1, 2 denote 
the species) 



= -7i (Vm - VtrQy(i)) + Cai{t) 



(54) 



where Vi is the particles velocity \'tray{t) = Aosin(a;t) 
is oscillating tray velocity, 7 provides linear damping, 
and Ccti is Gaussian white noise acting independently on 
each disk. The model reproduced segregation instabil- 
ity and subsequent coarsening of stripes. More realis- 
tic discret£_ekmentsnn^ were recently performed 
bv lCiamarra et al\ l)2005j) . In these simulations a binary 
mixture of round disks of identical sizes but two differ- 
ent frictions with the bottom plate (in fact, velocity- 
dependent viscous drag was assumed), separated in al- 
ternating bands perpendicular to the oscillation direc- 
tion irrespectively on initial conditions: both random 
mixed state and separated along the direction of oscil- 
lations state were used. Using particles of the same size 
eliminated the thermodynamic "excluded volume" mech- 
anism for segregation, and the authors argued that the 
mechanism at work is related to the dynamical shear in- 
stability similar to the Kelvin-Helmholtz instability in 
ordinary fluids. It was confirmed by a numerical obser- 
vation of the interfacial instability when two monolay- 
ers of grains with different friction constant were placed 
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in contact along a flat interface parallel to the direc- 
tion of horizontal oscillations. Similar in stability is ap- 
parently responsible for ripp le formation l)Scherer et all 
11999"; "Stcgncr and Wcsfrcid I. [T999|) . 

■Poolcy and Yeomans (,,200^ proposed theoretical de- 
scription of this experiment based on continuum model 
for periodically-driven two isothermal ideal gases which 
interact through frictional force. It was shown analyt- 
ically that segregated stripes form spontaneously above 
critical forcing amplitude. While the model reproduces 
the segregation instability, apparently it does not exhibit 
coarsening of stripes observed in the experiment. More- 
over, applicability of the isothermal ideal gas model to 
this experiment where the particles are almost at rest is 
an open question. 

Similar coarsening effect in granular segregation 
in a particula r ly sim ple geometry was studied by 
lAumaitre et all l)200l|) . They investigated the dynamics 
of a monolayer of grains of two different sizes in a dish 
shaken in a horizontal "swirling" motion. They observed 
that large particles tend to aggregate near the center of 
the cavity surrounded by small particles. The qualitative 
explanation of this effect follows from simple thermody- 
namic considerations (see above). Indeed, direct tracing 
of particle motion showed that the pressure in the area 
near the large particles is smaller than outside. But small 
particles do not follow the gradient of pressure and as- 
semble near the center of the cavity because this gradient 
is counterbalanced by the force from large particles. The 
inverse of force acting on large particles l eads to their ag- 
gregat ion near the center of the cavity. lAumaitre et ol\ 
J2OOI') proposed a more quantitative model of segregation 
based on the kinetic gas theory and found satisfactory 
agreement with experimental data. 

iBurtallv et al\ (|2002*) studied spontaneous separation 
of vertically vibrated mixtures of particles of similar sizes 
but different densities (bronze and glass spheres). At 
low frequencies and at sufficient vibrational amplitudes, 
a sharp boundary between the lower layer of glass beads 
and the upper layer of the heavier bronze spheres was 
observed. At higher frequencies, the bronze particles 
emerge as a middle layer separating upper and lower glass 
bead layers. The authors argue that the effect of air on 
the granular motion is a relevant mechanism of particle 
se paration. A somewha t similar conclusion was achieved 
by iMobius et al\ (^2001*) in experiments with vertically- 
vibrated column of grains containing a large "intruder" 
particle. 

lArndt et oi.l ^05); 'Fi odor and Ottinol 12003") per- 
formed detailed experiments on axial segregation in slur- 
ries, or bi-disperse grain-water mixtures. A mixture of 
two types of spherical glass beads of two sizes were placed 
in a water- filled tube at the volume ratio 1:2. Authors 
found that both rotation rate and filling fraction play an 
important role in band formation. Namely, bands are 
less likely to form at lower fill levels (20-30%) and slower 
rotation rates (5-10 rpm). They mostly appear near the 
ends of the drum. At higher fill levels and rotation rates. 



bands form fast er, and there are mo r e of them through- 
out th e drum. lArndt et aJ] ^QQ^; iFiodor and Ottinol 
linnS) ci'lso studied the relation between the bands visi- 
ble on the surface, and the core of small beads, and found 
that for certain fill levels and rotation speeds, the core re- 
mains prominent at all times, while in other cases the core 
disappears completely between bands of small particles. 
They also observed an interesting oscillatory instability 
of interfaces between bands at high rotation speeds. All 
these phenomena still await theoretical modelling. 



VIII. GRANULAR MATERIALS WITH COMPLEX 
INTERACTIONS 

A. Patterns in solid-fluid mixtures 

Presence of interstitial fluid significantly complicates 
the dynamics of granular materials. Hydrodynamic flows 
lead to the viscous drag and anisotropic long-range inter- 
action between particles. Even small amounts of liquid 
leads to cohesion among the particles which can have 
a profound effect on macroscopic properties of granu- 
lar assemblies such as angle of repose, avalanching, abil- 
ity to segregate, e tc. (see for example Sec. IVII Bland 
Li and McCarthvl ll2 005): S amadani and KiidrollJ ^20001 
200lj) : lTegzes et a^J ((2002,) ). 

In this Section we will discuss the case when the vol- 
ume fraction of fluid in the two-phase system is large, 
and the grains are completely immersed in fluid. This 
is relevant for many industrial applications, as well as 
for geophysical problems such as sedimentation, erosion, 
dune migration, etc. 

One of the most technologically important examples 
of particle-laden flows is a fluidized bed. Fluidized beds 
have been widely used since German engineer Fritz Win- 
kler invented the first fluidized bed for coal gasifica- 
tion in 1921. Typically, a vertical column containing 
granular matter is energized by a flow of gas or liquid. 
Fluidization occurs when the drag force exerted by the 
fluid on the granulate exceeds gravity. A uniform flu- 
idization, the most desirable regime for most industrial 
applications, turns out to be prone to bubbling insta- 
bility: bubbles of clear fluid are created spontaneously 
at the bottom, tra verse the granu lar layer and destroy 
the uniform state l|JacksonL |2000|) . Instabilities in flu- 
idized beds is an activ e area of research in the engi- 
neering community, se e llOida,sr)Owl [T99i I.Ta,cksonL 120001 
iKunii and L evenspiel Il99l|) . A shallow fluidized bed 
shows many similarities with mechanically vibrated lay- 
ers, see Se ction [Vl In par ticular, modulations of airflow 
studied bv lLi et oil l)2003|) result in formation of subhar- 
monic square and stripe patterns (see Fig. I23II similar 
to those in mechanically- vibra t ed sy stems (fMelo et all . 
ll994Lll99,'itfUmbanhowar et all\im(t\ . 

Wind and water driven granular flows play important 
roles in geophysical processes. Wind-blown sand forms 
dunes and beaches. The first systematic study of air- 
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borne (or aeolian) sand transport was conducted by R. 
Bagnold during Wold War II, see iBagnoldl 1119541) who 
identified two primary mechanisms of sand transport: 
saltation and creep, and proposed the first empiric re- 
lation for the sand flux q driyen by wind shear stress t: 



q = C, 




(55) 



where Cb = const ^ Ua is air density, d is the grain di- 
ameter, D = 0.25 mm is a reference grain size, and 
M* = yr/va is wind friction velocity. La ter many refine- 
ments of Ea. (|^ were proposed, see e.g. (|Pve and Tsoail 

11991^. 

iNishimori and Ouchil l)l993|) proposed a simple theory 
which describes formation of ripples as well as dunes. 
The theory is based on a lattice model which incorpo- 
rates separately saltation and creep processes. The model 
operates with the height of sand at each lattice side at 
discrete time n, hn(x, y). The full time step includes two 
substeps. The saltation substep is described as 



hn{x,y) 
hn{x + L{h{x,y)),y) 



hnix,y)~q (56) 
hn(x + L{hn{x,y)),y) +q 



where q is the height of grains being transferred from 
one (coarse grained) position (z, y) to the other position 
{x-\- L,y) on the lee side (wind is assumed blowing in the 
positive X direction) , L is the flight length in one saltation 
which characterizes the wind strength. It is assumed that 
q is conserved. Since the saltation length L depends on 
multiple factors, the following simple approximation is 
accepted 



L = Lo + bhnix,y) 



(57) 



with Lq measuring wind velocity and b = coiist. The 
creep substep involves spatial averaging over neighboring 
sites in order to describe the surface relaxation due to 
gravity. 



hn+i{x,y) = hn{x,y) + 



(58) 
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where X^jVAf ^"^^ TIiNNN denote summation over the 
nearest neighbors and next nearest neighbors corre- 
spondingly, and D — const is the surface relaxation 
rate. Despite its simplicity, simulation of the model 
reproduced formation of ripples and consequently ar- 
rays of barchan (crescent shaped) dunes, see Fig. ^2 
[Nishimori and Ouchi ( 1993) found that above certain 
threshold an almost linear relation holds between the 
selected wavelength of the dune pattern and the "wind 
strength" L. 

In the long- wave limit Eqs. H57|) - (|59|l can be reduced 
to more traditional continuum models considered below. 



FIG. 41 Sand ripple pattern (top panel) and barchan dunes 
(lower panel) obtained from simulations of Eqs. 15711 - 159(1 
INishimori and Ouchil lll993l) 



FIG. 42 Inte raction and coars ening of one-dimensional dune 
system, from iPrigozhinI lll999f) 



In the continuum description of the evolution of the 
sand surface, the profile h is governed by the mass con- 
servation equation 



Vsdth = -Vq, 



(59) 



where is the density of sand and q is the sand flux. 
In order to close Eqs. (|S51) . lf5^ . several authors pro- 
posed different phenomenological relations between shear 
stress at the bed surface r and the height h, see e.g . 
llAndreotti et^ali \2m 2: Hersen et ai, 2004; Krov et all 
l2002a.k .Nishimori and Ouchiill993tiPrigozhinlll999() . 

There are ma ny theories generalizing 

iNMnmoryind Ouchi (1993*) app roach, se e e.g . 
l)CaDS anTvaanderwalle. ,2001^1 . jPrigozhinl l)l999|) 
described the evolution of dunes by a system of two 
equations simi lar to the BC RE model discussed earlier 
in Sec. EH (,Bouchaud et all J994, 1995) . One equation 
describes the evolution of the local height h while 
another equation describes the density R of particles 
rolling above the stationary sand bed profile (reptating 
particles) , 



dth = T{h, R)-f 

dtR = -VJ + Q - r(/i, i?) 



(60) 
(61) 



where F is the rolling-to-steady sand transition rate, J 
is the horizontal projection of the flux of rolling parti- 
cles, Q accounts for the influx of falling reptating grains, 
and / is the erosion rate. With an appropriate choice 
of rate functions F, /, Q and J, Eqs. H61|) can reproduce 
many observed features of dune formation, such as initial 
instability of flat state, asymmetry of the dune profiles, 
coarsening and interaction of dunes, etc., see Fig . 

Thus, simplified n i odels such as llKrov et alV l2002at 
INishimori and Ouchil [l993t jPrigozhinl Il999|) have been 
successful in explaining many features of individual dune 
growth and evolution, see Fig. 021 However we should 
note that up to date none of the dune models have been 
able to address satisfa ctorily the waveleng th selection in 
large-scale dune fields ijHersen et al\ . l2004j) . 

The phenomenon qualitatively similar to the dune for- 
mation occurs in an oscillatory fluid flow above a gran- 
ular layer: sufficiently strong flow oscillations produce 



FIG. 43 Evolution of two barchan dunes described by Eqs. 
15511 . 1591 1. Small dune i s und ersupplied and eventually 
shrinks, from .Hersen et all i2004l) 
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FIG. 44 Fingering instability of planar avalanche in the un- 
derwater flow, figure on the left zooms on individual finger, 
from .Malloggi et ali ^2005a^ 



so-called vortex ripples on the surface of the underlying 
granular layer. These ripples are familiar to an y beach- 
goer. V ortex ri pple formation was first stud ied bvlAvrtor 
1 191C Y. iBaenol d (195 (1), and recent ly bv IScherer et al 
l l999D :ISteg ner and Wesfreid ] (119991) ' and others. It was 
found that ripples emerge via a hysteretic transition, and 
are characterized by a near-triangular shape with slope 
angles close to the repose angle. The characteristic size 
of the ripples A is directly proportional to the displace- 
ment amplitude of the fluid flow a (with a proportionality 
constant w 1.3) and is roughly independent on the fre- 
q uency. 

lAndersen et all (|200lf) introduced order parameter 
models for describing the dynamics of sand ripple pat- 
terns under oscillatory flow based on the phcnomcnolog- 
ical mass transport law between adjacent ripples. The 
models predict the existence of a stable band of wave 
numbers limited by secondary instabilities and coarsen- 
ing of small ripples, in agreement with experimental ob- 
se rvations. 

iLanglois and Valancd (|2005() studied underwater rip- 
ple formation on a two-dimensional sand bed sheared by 
viscous fluid. The sand transport is described by gener- 
alization of Eq. (|55|) taking into account both the local 
bed shear stress and the local bed slope. Linear stabil- 
ity analysis revealed that ripple formation is attributed 
to a growing longitudinal mode. The weakly nonlin- 
ear analysis taking into account resonance interaction of 
only three unstable modes revealed a variety of steady 
two-dimensional ripple patterns drifting along the flow 
at some speed. 

Experiments i n dune formation h ave been recently per- 
formed in water ("Beta t et ai\.\l99^ . While water-driven 
and wind driven dunes and ripples have similar shape, 
the underlying physical processes are likely not the same 
due to a different balance between gravity and viscous 
drag in air and water. 

Spectacular erosion patterns in sediment granular lay- 
ers were observed in ex periments with underwater flows 
(|Daerr et all l2003t l^aUoggi et at, 20053). In particu- 
lar, a fingering instability of flat avalanche fronts was 
observed, see Fig. ^] These patterns are remarkably 
similar to those in thin films on inc lined surfaces , both 
with clear and pa rticle-laden fluids (jTroian et all Il989t 
IZhou et adl2005f) . In the framework of lubrication ap- 
proximation the evolution of fluid film thickness h is de- 
scribed by the following dimensionless equation following 
from the mass conservation law: 



dth + V ■{ [/i^VV^/i] - Dh^Vh] + d^h^ = (62) 

where dimensionless parameter D is inversely propor- 
tional to water surface tension. The instability occurs 



FIG. 45 Upper panel: Self-supporting knolls formed in wa- 
ter/glass beads suspension in a horizontally rotating cylinder 
(side and end views). Lower panel: Computational results: 
schematics of the flow (a); the height of computed knoll struc- 
tures (b,c), from iDuong et al\ J2004I) 



for small D values, i.e. in the large surface tension limit. 
However, despite visual similarity the physical mecha- 
nism leading to this fingering instability is not obvious: 
in fiuid films the instability is driven (and stabilized) by 
the surface tension, whereas in the underwater granular 
fl ow fluid surfa c e tens ion plays no role. 

iDuong et all ()2004|) studied formation of periodic ar- 
rays of knolls in a slowly rotating horizontal cylinder filled 
with granular suspension, see Fig. 031 The solidified sed- 
iment knolls co-exist with freely circulating fluid. The 
authors applied variable viscosity fluid which formally 
allows simultaneous treatment of solid and liquid phase. 
In this model the effective flow viscosity /Xs diverges at 
the solid packing fraction 4'rcp, 
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(63) 



where /io is the clear fiuid viscosity and b is an empir- 
ical coefficient. The model qualitatively reproduces the 
experiment, see Fig. 1451 An interesting question in this 
context is whether there is a connection to the experi- 
ment by Shcn ( 2002) where somewhat similar structures 
were obtained for the flow of "dry" particles in a hori- 
zontally rotating cylinder. 

As it was mentioned in Sec. IVICT IConwav et a,l\ 
()2004|) reported that an air-fluidized vertical column 
of bi-disperse granular media sheared between counter- 
rotating cylinders exhibits formation of nontrivial vortex 
structure strongly reminiscent o f Taylor vortices in con- 
ventional fluid, see e.g. (Andcrcc k et al 1IT986). Authors 
argue that vortices in fluidized granular media, unlike 
Taylor vortices in fluid, are accompanied by the novel 
segregation-mixing mechanism specific for granular sys- 
tems, see Fig. QU] Interestingly, no vortices were ob- 
served in a similar experi ment in Couet t e geom etry with 
m onodisperse gl a ss bea ds (|Losert et a/.l . l2000() . 

Ilvanova et all (Il996() studied patterns in a horizontal 
cylinder filled with sand/liquid mixture and subject to 
horizontal vibration. For certain vibration parameters 
standing wave patterns were observed at the sand/liquid 
interface. Authors argue that these wave patterns are 
similar to the Faraday ripples found at liquid/liquid in- 
terface under vertical vibration. 



B. Vortices in vibrated rods 

In Section we reviewed instabilities and collective 
motion in mechanically vibrated layers. In most exper- 
iments the particle shape was not important. However, 
strong particles anisotropy may give rise to non-trivial 
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FIG. 46 Phase diagram for the system of verticaUy-vibrated 
rods, driving frequency 50 Hz. Vortices are observed for suffi- 
ciently high fiUing f r action n/ and above critical acceleration 
r, from lBlair et al\ ll2003al) 



FIG. 47 Azimuthal velocity of the vortex vs distance from the 
center for different parameter values, from iBlair et aL i.2003al 

effects. IVillarruel et all l|200Cfl observed onset of nematic 
order in packing of long rods in a narrow vertical tube 
subjected to vertical tapping. The rods initially com- 
pactify into a disordered state with predominantly hori- 
zontal orientation, but at later times (after thousands of 
taps) they align vertically, first along the walls, and then 
throughout the volume of the pipe. The nematic order- 
ing can be understood in terms of the excluded volume 
a rgument put forwar d by Dnsa ger (1949). 

iBlair et aj} |2003a|l studied the dynamics of vibrated 
rods in a shallow large aspect ratio system. Surprisingly, 
they found that vertical alignment of rods at large enough 
filling fraction n f and the amplitude of vertical accelera- 
tion (r > 2.2) can occur in the bulk, and it does not re- 
quire side walls. Eventually, most of the rods align them- 
selves vertically in a monolayer synchronously jumping 
on the plate, and engage in a correlated horizontal mo- 
tion in the form of propagating domains of tilted rods, 
multiple rotating vortices etc, see Fig. E| and Fig. HSl 
The vortices exhibit almost rigid body rotation near the 
core, and then the azimuthal velocity falls off. Fig. W7\ 
The vortices merge in the course of their motion, and 
eventually a single vortex is formed in the cell. 

Experiments showed that the rod motion occurs when 
the rods are tilted from the vertical, and it always 
occurs in the direc tion of tilt. In subsequent work 
IVolfson et al\ l|2004f) experimentally demonstrated that 
the correlated transport of bouncing rods is also found 
in quasi-one-dimensional geometry, and explained this 
effect using molecular dynamics simulations and a de- 
tailed description of inelastic frictional contacts be- 
tween the rods and the vibrated plate. Effectively, 
bouncing rods become self-propelled objects similar 
to other self-propelled systems, for which large-scale 
coherent motion is often observed (bird flocks, fish 
scho ols, chemotactic microorgani s m aggregation, e tc., se e 
e.g. iGregoire and Chati ll2004l): iHelbing et al\ l|2000|) : 
'Helbing ('2001'); 'Toner and Tu' ('1995*)). 

Aranson and Tsimring ( 2003) developed a phe- 
nomcnological continuum theory describing coarsening 
and vortex formation in the ensemble of interacting 
rods. Assuming that the motion of rods is overdamped 
due to the bottom friction, the local horizontal velocity 
V = {vxTVy) of rods is of the form 

^ = -{Vp-anh{n)v) /Cv, (64) 

where v is the density, p is the hydrodynamic pressure, 
the tilt vector n = {nx,ny) is the projection of the rod 



FIG. 48 Sequence of images illustrating coalescence of vor- 
tices in the model of vibrated rods, the field |n| is shown, 
black dots corresponds to vo rtex cores where |n| = 0, from 
lAranson and Tsimrind (120031) 

director on the {x, y) plane normalized by the rod length, 
i.e n = |n |, and C is friction coe f ficient . According to 
iBlair et al\ l|2003aj) : IVolfson et all l|2004|) . the rods drift 
is determined by the average tilt of neighboring rods, 
thus the term an.fQ{n)v accounts for the average driving 
force from the vibrating bottom on the tilted rod. Eq. 
(|64|l combined with the mass conservation law yields 

dtv = — div(v!^) = C "'^div (Vp — anfQ{n)v) . (65) 



To account for the experimentally observed phase t 
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review) by assuming that pressure p can be obtained from 
the variation of a generic bistable "free energy" func- 
tional F with respect to the density field p = 5F/5v. 

To close the description the equation for the evolution 
of tilt n is added on generic symmetry arguments: 

dtn = /i(i.)n- Inpn-H 

+ f2{v) (CiV^n + ^Vdivn) + p^v. (66) 

Here /i^2 are certain functions of v, ^1^2 characterize dif- 
fusion coupling between the neighboring rods. Since the 
tilt field is not divergence-free, from the general symme- 
try considerations both ^1.2 7^ ^. 

Numerical and analytic studies of Eqs. (|65() . (|66|l re- 
vealed phase coexistence, nucleation and coalescence of 
vortices in accord with the experiment, see Fig. HSl 

An interesting experi ment with anisotropic chiral par- 
ticles was performed bv iTsai ~ a/. (2005). The role of 
particles was played by bend-wire objects which rotated 
in a preferred direction under vertical vibration. The ex- 
periments demonstrated that individual angular rotation 
of the particles was converted into the collective angular 
momentum of the granular gas of these chiral objects. 
The theoretical description of this system was formulated 
in the framework of two phenomenological equations for 
the density v and center-of-mass momentum density I'v 
and the spin angular momentum density I — Xfl aris- 
ing from the ration of particles around their center of 
mass, is the particle's rotation frequency. Whereas the 
equations for density and velocity are somewhat similar 
to those for the vibrated rod system, the equation for 
the spin momentum clearly has no counterpart in the vi- 
brated rod system and was postulated in the following 
form: 

dtl + Vwl=T- F^^ - r{n -uj) + DnV^n (67) 



^ These constants are analo gous to the first and second viscosity 
in ordinary fluids, see e.g. iLandau and Lifshitst^l959^ 
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where t is the source of the angular rotation (due to chi- 
rahty of particles), uj is coarse-grained or collective an- 
gular velocity, and F are dissipative coefficients due 
to friction and Dji is the angular momentum diffusion. 
Eq. H67() predicts, in agreement with the experiment, the 
onset of collective rotation of the gas of particles. Possi- 
bly, it also exhibits non-trivial spatio-temporal dynamics 
similar to those in the system of vibrated rods. How- 
ever, due to the small number of particles (about 350) in 
the experiment the nontrivial collective regimes were not 
reported. 



C. Electrostatically driven granular media 

Large ensembles of small particles display fascinating 
collective behavior when they acquire an electric charge 
and respond to competing long-range electromagnetic 
and short-range contact forces. Many industrial tech- 
nologies face the challenge of assembling and separating 
such single- or multi-component micro and nano-size en- 
sembles. Traditional methods, such as mechanical vibra- 
tion and shear, are infective for very fine powders due to 
agglomeration, charging, etc. Electrostatic effects often 
change statistical prope rties of granular matter such as 
energy dissipation rate ijSheffier and Wolil2002l). veloc- 
ity distributions in granular gases ranson and OlafsenL 
|2P02; Kohlstcdt et ai, 2005), agglomeration rates in sus- 
p ensions dDamm er and Wolil2004) . etc. 

lAranson et al\ l|2000l l2002() : ISapozhnikov et al\ 
ll2003ai l200l " studied electrostatically driven granular 
matter. This method relies on the collective interactions 
between particles due to a competition between short 
range collisions and long-range electromagnetic forces. 
Direct electrostatic excitation of small particles offers 
unique new opportunities compared to traditional 
techniques of mechanical excitation. It enables one to 
deal with extremely fine nonmagnetic and magnetic 
powders which are not easily controlled by other means. 

In most experimental realizations, several grams of 
mono-dispersed conducting micro-particles were placed 
into a 1.5 mm gap between two horizontal 30 x 30 
cm^ glass plates covered by transparent conducting lay- 
ers of indium tin-dioxide. Typically 45 /im Copper 
or 120 /zm Bronze spheres were used. Experiments 
were also performed with much smaller 1 fim particles, 
l|SaDozhnikov et qiJ . l200l . An electric field perpendicu- 
lar to the plates was created by a high voltage source (0-3 
kV) connected to the inner surface of each plate. Experi- 
ments were performed in air, vacuum, or in the cell filled 
with non-polar weakly-conducting liquid. 

The basic principle of the electro-cell operation is as 
follows. A particle acquires an electric charge when it 
is in contact with the bottom conducting plate. It then 
experiences a force from the electric field between the 
plates. If the upward force induced by the electric field 
exceeds gravity, the particle travels to the upper plate, 
reverses charge upon contact, and is repelled down to the 



bottom plate. This process repeats in a cyclical fashion. 
In an air-filled or evacuated cell, the particle remains im- 
mobile at the bottom plate if the electric field E is smaller 
than the first critical field Ei. For E > Ei a.n isolated 
particle leaves the plate and starts to bounce. However, 
if several particles are in contact on the plate, screening 
of the electric field reduces the force on individual par- 
ticles, and they remain immobile. A simple calculation 
shows that for the same value of the applied electric field 
the force acting on isolated particles exceeds by a factor 
of two the force acting on the particle inside the dense 
monolayer. However, if the field is larger than a second 
critical field value, E2 > Ei, all particles leave the plate, 
and the system of particles transforms into an uniform 
gas-like phase. When the field is decreased below E2 
{El < E < E2), in air- filled or evacuated cells localized 
clusters of immobile particles spontaneously nucleate to 
form a static clusters (precipitate) on the bottom plate 
( Aranson et a i. 200(j). The clusters exhib it the Ostwald- 
type ripening ijBravl [l 993 iMeersonL 1 1 99(^ . see also Sub- 
sec, nvm 

1. Coarsening of clusters 

Results for the electrostatically driven system yielded 
the following asymptotic scaling law, see Fig. EOI 

N^^ (68) 

where N is the number of clusters and t is time. Accord- 
ingly, the average cluster area (A) increases with time as 
(A) ~ t. This behavior is c onsistent with t he interface- 
controlled Ostwald ripening llMeersonL Il99(j) . 

A theoretical description of coarsening in an elec- 
trostatically d r iven g r anular system was develo ped by 
lAranson et d] ()2000|) . ISapozhnikov et al\ ()2003|) . The 
theory was formulated in terms of the Ginzburg-Landau- 
type equation for the number density of immobile parti- 
cles (precipitate or solid) n 

dtn = \'^n + (t>{n,ng) (69) 

where ng is the number density of bouncing parti- 
cles (gas) Ug, and 4>{n,ng) is a function characterizing 
a solid/gas conversion rate. The effectiveness of the 
solid/gas transitions is controlled by the local gas con- 
centration Ug. It was assumed that the gas concentration 
is almost constant because the particle's mean free pass 
in the gas state is very large. The gas concentration Ug 
is coupled to n due to total mass conservation constraint 

Sug + j j n{x,y)dxdy — AI, (70) 

where 5* is the area of domain of integration, and M 
is the total number of particles. Function (j)(n,ng) is 
chosen in such a way as to provide bistable local dynamics 
of concentration corresponding to the hysteresis of the 
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FIG. 49 Illustration of phase separation and coarsening dy- 
namics. (a)-(c) Numerical solution of Eqs. 1691 . 17UII . white 
corresponds to dense cluste rs, black to dilute gas , (d)-(f) show 
experimental results, from lAranson et al\ 1I2OO2D . 



FIG. 50 Average cluster area {A{t)) (a) and inverse num- 
ber of clusters 1/N{t) vs time in air-filled cell. The straight 
line in (b) shows theor etical prediction 1/N ~ t, from 
ISapozhnikov et al\ ll2005t) 



gas/solid transition. The above description yields a very 
similar temporal evolution of clusters (see Fig. I49|l and 
produces a correct scaling for the number of clusters Eq. 
(EHIl- 

In the so-called sharp interface limit when the size of 
clusters is much larger than the width of interfaces be- 
tween clusters and granular gas, Eq. can be reduced 
to equations for the cluster radii Ri (assuming that clus- 
ters have circular form): 



dR, 
~dt' 



Rc{t) 



1 

Ri 



(71) 



where Rc is critical cluster size, k is effective surface ten- 
sion (experimental mea surements of cluste r surface ten- 
sion were conducted bv ISapozhnikov et al\ (|^003)). The 
critical radius Rc is a certain function of the granular gas 
concentration Ug that enters Eqs. (|71|l through the the 
conservation law Eq. (|70|l which in two dimensions reads 



N 

1=1 



(72) 



The statistical properties of Ostwald ripening can 
be understood in terms of the probability distri- 

bution function f(R ,t) of clust e r sizes. Follo wing 

I Lifshitz and Slvozovl ill 9581 llflfilh : IWagneil lll96lh and 

neglecting cluster merger, one obtains in the limit N — > 
00 that the probability distribution f{R,t) satisfies the 
continuity equation 



dtf + Or (i?/) = 0. 



(73) 



From the mass conservation in the limit of small gas con- 
centration Eq. (|72|l one obtains an additional constraint; 



/"OO 

TT / R^f{R,t)dR = M 
Jo 



(74) 



Eqs. 1)731) . 174|l have a self-similar solution in the form 



R 



(75) 



For the total number of clusters N = fdR the scal- 
ing Ea. (|75|l yields N ~ 1/t, which appears to be in a 
good agreement with the experiment, see Fig. 1501 How- 
ever, the cluster size distribution function appears to be 
in a strong disagreement, see Fig . ICT In parti cular, 
I Lifshitz and Slvozovl lll958l Il96ll) : IWagneil lll96lD the- 
ory predicts the distribution with a cut-off (dotted Hue) 
whereas the experiment yields the function with an ex- 
ponential tail. A much better agreement with the exper- 
iment was obtained when binary coalescence of clusters 



FIG. 51 Scaled cluster size distribution function F{^) with 
5 = R/Vt. The squares show experimental results, the dot- 
ted line shows a n alytic result form Lifshitz-Slyozov- Wagner 
theory iWagneil Il96ll) . and solid line shows F obtained 
from the theory accoun ting for binary coalescence, from 
ISapozhnikov et al\ ll2005l) 



was incorpor ated in the Lifshitz-Slyozov- Wa gner theory 
l|Conti et ad l2002: Sap ozhnikov et a/J . l2005j) .The coales- 
cence events become important for a finite area fraction 

o f the clusters. 

iBen-Naim and Krapivskvl l)2003|) applied an exchange 
growth model to describe coarsening in granular media. 
In this theory the cluster growth rates are controlled only 
by the cluster area ignoring shape effects. Assuming that 
the number of particles in a cluster evolves via uncorre- 
lated exchange of single particles with an other cluster 
the following equation for the density of clusters contain- 
ing k particles can be derived: 

^ = ^ AAjK,, {Sk,^+l + - 25k,^) (76) 
hi 

where A^. is the probability to find a cluster containing 
k particles, Kij the exchange kernel and dk,i is the Kro- 
necker symbol. For the choice of homogeneous kernel 
Kij = {ij)^ with A = 1 this theory predicts correct scal- 
ing of the cluster size with time R ^ ^/t and exponential 
decay of the cluster size distribution function, as in the 
experiment. The choice of A = 1 is equivalent to the 
assumption that the exchange rate is determined by the 
size o f the cluster. In the theory by ISapozhnikov et a/l 
lioO^ the cluster evolution is governed by the evapora- 
tion/deposition of particles at the interface of the cluster 
and controlled by the overall pressure of the granular gas. 
Thus, both theories predict the same scaling behavior, 
however the underlying assumptions are very different. 
A possible explanation for this may be that while the ex- 
change growth model ignores the curvature of the cluster 
interface and the dependence on exchange rate on the 
pressure of granular gas, the agreement is obtained by 
tuning the adjustable parameter A. 



2. Dynamics of patterns in a fluid-filled cell 



ISapozhnikov et all l)20n3a() performed experiments 
with electrostatically driven granular media immersed 
in a weakly conducting non-polar fluid (toluene-ethanol 
mixture). Depending on the applied electric field and 
the ethanol concentration (which controls the conductiv- 
ity of the fluid), a plethora of static and dynamic patterns 
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were discovered, see Fig. For relatively low concen- 
trations of ethanol (below 3%), the qualitative behavior 
of the liquid-filled cell is not very different from that of 
the air-filled cell: clustering of immobile particles and 
coarsening were observed between two critical field val- 
ues Ei^2 with the clusters being qualitatively similar to 
that of the air cell. However, when the ethanol concen- 
tration is increased, the phase diagram becomes asym- 
metric with respect to the direction of the electric field. 
Critical field magnitudes, E12, are larger when the elec- 
tric field is directed downward on the upper plate) 
and smaller when the field is directed upward ("— " on 
the upper plate). This difference increases with ethanol 
concentration. The observed asymmetry of the critical 
fields is apparently due to an excess negative charge in 
the bulk of the liquid. 

The situation changes dramatically for higher ethanol 
concentrations: increasing the applied voltage leads to 
the formation of two new immobile phases: honeycomb 
(Fig. If 5b l for the downward direction of the applied 
electric field, and two-dimensional crystal-type states for 
the upward direction. 

A further increase of ethanol concentration leads to the 
appearance of a novel dynamic phase - condensate (Fig. 
If 5b .d) where almost all particles are engaged in a circular 
vortex motion in the vertical plane, resembling Rayleigh- 
Benard convection. The condensate co-exists with the di- 
lute granular gas. The direction of rotation is determined 
by the polarity of the applied voltage: particles stream 
towards the center of the condensate near the top plate 
for the upward field direction and vice versa. The evo- 
lution of the condensate depends on the electric field di- 
rection. For the downward field, large structures become 
unstable due to the spontaneous formation of voids (Fig. 
If 51 FI. These voids exhibit complex intermittent dynam- 
ics. In contrast, for the upward field, large vortices merge 
into one, forming an asymmetric object which often per- 
forms composite rotation in the horizontal plane. The 
pattern formation in this system is most likely caused 
by self-induced electro-hydrodynamic micro- vortices cre- 
ated by the particles in weakly-conducting fluids. These 
micro-vortices create long-range hydrodynamic vortex 
flows which often overwhelm electrostatic repulsion be- 
tween likely-charged particles and introduce attractive 
dipole-like hydrodynamic interactions. Somewhat simi- 
lar micro- vortices are known in driven colloidal systems, 
se e e.g. (Ych et al, I997j. 

lAranson and Sapozhnikovl l)2004|) developed a phe- 
nomenological continuum theory of pattern formation for 
metallic micro-particles in a weakly conducting liquid 
subject to an electric field. Based on the analogy with the 
previously deve l oped theory of coarsening in air-field cell 
i)Aranson eF^ . 120021) . the model is formulated in terms 
of conservation laws for the number densities of immobile 
particles (precipitate) Up and bouncing particles (gas) Ug 
averaged over the thickness of the cell: 

dtup = VJp + / , dtug = VJg - /. (77) 



FIG. 52 Sequence of snapshots illustrating evolution of pul- 
sating rings (top raw) and rotating vortices (bottom raw) 
obtained from numerical solutio n of Eqs. 177L179j . from 
lAranson and Sapozhnikovl i2004) 

Here Jp^g are the mass fluxes of precipitate and gas re- 
spectively and the function / describes gas/precipitate 
conversion which depends on np,g, electric fleld E and 
local ionic concentration c. The fluxes are written as: 

Jp,3 = Dp^gVup^g + ap^g{E)vj_np^g{l - P{E)np^g), (78) 

where v± is horizontal hydrodynamic velocity, Dp^g are 
precipitate/gas diffusivities. The last term, describ- 
ing particles advection by fluid, is reminiscent of the 
Richardson-Zaki relation for a drag force frequently 
used in the engineering literature (Richardson and ZakJ, 
If954j) . The factor (I — P{E)np^g) describes the satura- 
tion of flux at large particle densities n ~ due to the 
decrease of void fraction. Terms ap,g describe advec- 
tion of particles by the fluid. Interestingly, in the limit of 
very large gas diffusion Dg ^ Dp and without advection 
terms (ofp.g = 0) the mod el reduces to Eq s. (1^ and l(7n|) 
applied for air-filled cell ijAranson et a/.l . l2b02i) . 

Eqs. (|77|l are coupled to the cross section averaged 
Navier-Stokes equation for vertical velocity Vz- 

no{dtv,+vWv,)= fiW^v,-d,p + Ezq (79) 

where hq is the density of liquid (we set hq = f), /i is 
the viscosity, p is the pressure, and q is the charge den- 
sity. The last term describes the electric force acting on 
charged liquid. Horizontal velocity v± is obtained from 
Vz using the incompressibility condition dzV^ + V _\_vj_ =^ 
in the approximation that vertical vorticity f2z ~ dxVy — 
dyVx is small compared to in-plane vorticity. This as- 
sumption allows one to find the horizontal velocity as a 
gradient of quasi-potential (p: = — Vj^0. 

For an appropriate choice of the parameters the model 
Eqs. ((77|l . (|7^ yields qualitatively correct phase diagram 
and the patterns observed in the experiment, see Figs. 
iniandEll 

D. Magnetic particles 

Electric and magnetic interactions allow introduc- 
tion o f controlled long-r a nge forces in granular sys- 
tems. iBlair et al\ ll2003ajl: iBfair and Kud roUi (2003b); 
IStambaugh et a L (2004aD) performed experimental 
studies with vibrofluidized magnetic particles. Several 
interesting phase transitions were reported, in particu- 
lar, the formation of dense two-dimensional clusters and 
loose qu asi-one-dimensional chains and rings. Blair et q3 
l)2003a|) considered pattern formation in a mixture of 
magnetic and non-magnetic (glass) particles of equal 
mass. The glass particles played the role of "phonons", 
their concentration allowed an adjustment of the typi- 
cal fluctuation velocity of the magnetic subsystem. The 
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FIG. 53 Phase diagram illustrating various regimes in mag- 
netic granular media, T is the temperature determined from 
the the width of velocity distri bution, <& is surfac e covera ge 
fraction of glass particles, from lBlair and KudroUil (l2003btl 

phase diagram delineating various regimes in this system 
is shown in Fig. While the phase diagram shows 

some similarity with equilibrium dipolar fluids (such as 
phase coexistence), most likely there are differences due 
t o the non-equilibrium char acter of granular systems. 

IStambaugh et all l)2004a|) performed experiments with 
relatively large particles (about 1.7 cm), and near the 
closed-packed density. It was found that particles form 
hexagonal closed-packed clusters in which the magnetic 
dipoles lay in the plane and assume circulating vorti- 
cal patterns. For lower density ring patterns were ob- 
served. Experiments with mixture of particles with two 
different magnetic moments revealed segregation effects 
jStambaugh et al, 20 04b|) . The authors argue that the 
static configurational magnetic energy is the primary fac- 
tor in pattern selection. 

Experiments bv i Blair and KudroUj l)2003b(l : 
IStambaugh et al\ l|2004albD were hmited to a small 
number (about 10"^) of large particles due to the intrinsic 
limitation of the mec hanical vibrofluidization technique. 
ISnezhko et al\ ()2005 l) performed experimental studies of 
90 iJ,m Nickel micro-particles subjected to electrostatic 
excitation, see also Subsec. IVIII.CI The electrostatic 
system allowed researchers to perform experiments with 
a very large number of particles (of the order of 10^) 
and a large aspect ratio of the experimental cell. Thus 
the transition between small chains and large networks 
(Fig. [T^ was addressed in detail. An abrupt divergence 
of the chain length was found when the frequency of 
field oscillations decreased, resulting in the formation of 
a giant interconnected network. 

Studies of the collective dynamics and pattern forma- 
tion of magnetic particles are still in the early phases. 
While it is natural to assume that magnetic interaction 
plays a dominant role in pattern selection, further com- 
putational and theoretical studies of pattern formation 
in systems of driven dipolar particles are necessary. Be- 
sides a direct relevance for the physics of granular media, 
studies of magnetic granular media may provide an ad- 
ditional insight into the behavior of dipolar hard sphere 
fluids wh ere the nature of solid/liquid transitions i s still 
debated l)de Gennes and Pincud. llQTOt iLevinL Il999|) . Vi- 
bration or electrostatically fluidized magnetic particles 
can also be viewed as a macroscopic model of a ferrofluid, 
where similar experiments are technically difhcult to per- 
form. 



IX. OVERVIEW AND PERSPECTIVES 

Studies of granular materials are intrinsically inter- 
disciplinary and they borrow ideas and methods from 



other fields of physics such as statistical physics, me- 
chanics, fluid dynamics, and the theory of plasticity. On 
the flip side, progress in understanding granular mat- 
ter can be often applied to seemingly unrelated phys- 
ical systems, such as ultra-thin liquid films, foams, 
colloids, emulsions, suspensions, and other soft con- 
densed matter systems. The common feature shared 
by these systems is the discrete microstructure di- 
rectly influencing macroscopic behavior. For example, 
the order parameter description similar to that of Sec. 



IVi.A.l 


was 


applied to stick-slip friction in ultra-thin 


films, ( 


Aranson et at. 


2002c; 


Carlson and Batistalll996l 


llsraelachvih et aZ.'. '198& .Urbach et alU2004}. 


tLomaitro 


(2002): iLemaitre and Carlson 1 l|2004|) ao- 



plied the idea of shear-transform ation zone (STZ) pio- 
neered by iFalk and Langeil l)l998|) for amorphous solids 
both to granular matter and to the boundary lubrica- 
tion problem in confined fluid. In this theory the plas- 
tic deformation is represented by a population of meso- 
scopic regions which may undergo non-affine deforma- 
tions in response to stress. Concentration of STZs in 
amorphous material is somewhat similar to the order pa- 
rameter (relative conc entration of defects) introduced by 
lAranson et al\ ll2002cll. A co nceptually similar approach 
was proposed bv lStaron et a l. (2002) who described the 
onset of fluidization as a percolation of the contact net- 
work with fully mobilized friction. Whereas derivation of 
the constitutive relations from first-principle microscopic 
rules is still a formidable challenge, these approaches are 
promising for understanding of not only the boundary 
lubrication problem, but also onset of motion in dense 
granular matter. 

Flowing liquid foams and emulsions share many simi- 
larities with granular matter: they have internal discrete 
structure (bubbles and drops play the role of grains) , and 
two different mechanisms are responsible for the trans- 
mission of stresses: elastic for small stress and visco- 
plastic above certain yield stress. However, there are 
additional complications: bubbles are highly deformable 
and, unlike granular matter, a number of particles may 
change due to the coalescence of bubbles. 

Foams and granular materials often exhibit sim- 
ilar behavior, such as non-trivial stress relaxation 
and power-law distr ibution of rearrangement events 
(Dcnnin and K nobleiill9 97). Sti ck-slip behavior was re- 
ported both for sheare d foams i Laurids en et all 120021) 
and granular materials l)Nasuno aU il997j) . Remark- 
ably, recent experimen ts with two dimensional foams 
( Lauridsei^et aZ., 2004) a nd three d imensional emulsions 
( Coussot et al . 20J12a b: iDaCruz I 12002.) strongly sug- 
gest the coexistence between flowing (liquid) and jammed 
(solid) states reminiscent of that in granular matter. 
Furthermore, avalanche behav ior reminiscent of granular 
flows down an i nclined plane llDaerr and Douadvl Il999|) 
was reported bv lCoussot et a.l\ ll2002al) for clay suspen- 
sions, see Fig. There are many approaches treat- 
ing foams, gel and suspensions as complex fluids with 
specific stress-strain constitutive relation. For example. 
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FIG. 54 Sequence of snapshots illustrating ev olution of clay 
suspen sion drop poured over sandpaper, from ICoussot et al\ 
^2002a^ 

FIG. 55 A possible phase diagram for jamming. The jammed 
region, near the origin, is enclosed by the depicted surface. 
The line in the temperature-load plane is speculative, and 
indicates how the yield stress might va ry for jammed systems 
in which there is thermal motion, from lLiu and Naeell J1998H 

iFuchs and Gated l)2002j) used the analogy between glasses 
and dense colloidal suspensions and applied the mode 
coupling approach to understand the nonlinear rheology 
and yielding. Similar approaches can be possibly usefu l 
fo r granular mate r ials llSchofield and OppenheimL Il994() . 

iLiu and Naeell l)l998|) suggested that a broad class 
of athermal soft matter systems (glasses, suspensions, 
granular materials) shows a universal critical behavior 
in the vicinity of solid-fluid or jamming transition, see 
Fig. |55l Whether jammed systems indeed have com- 
mon features that can be described by a universal phase 
diagram is an open issue. An interesting question in 
this context is a possibility of thermodynamic descrip- 
tion of driven, macroscopic, athermal systems like gran- 
ular materials and foams in terms of some kind of effec- 
tive te mperature. Studies of interacting particles under 
shear dCorvin a/.l . l2005t iMakse and Kurchanl 12002 : 
1 'Hern et all 120041: lOno et all 120021: IXu and O'HerE . 
12005") indicate that indeed under certain conditions it is 
possible to define an effective temperature (for example, 
from the equivalent of the Einstein-Stokes relation) for a 
broad class of athermal systems from comparison of the 
mechanical linear response with the corresponding time- 
dependent fluctuation-dissipation relation. However, the 
possibility of developing nonequilibrium thermodynamics 
of the basis of the effective temperature is under debate. 

Granular systems exhibit many similarities with traf- 
fic fiows and collective motion of self-propelled particles 
such as swimmi ng bacte r ia, fis h schools, bird flocks, etc., 
see for review (iHelbingl l2(ffll . In particular, jamming 
transition in granular media and traffic jams show simi- 
lar features, such as hysteresis, and clusters formation. 
Moreover, continuum models of traffic fiows are often 
cast in the form of modified Navier-Stokes equation with 
density-dependent viscosity, similar to granular hydrody- 
namics. 

Let us discuss briefly some open questions in the 
physics of granular matter. 

• Static vs. dynamic description. Gommonly ac- 
cepted models of rapid granular flows (granular hy- 
drodynamics) and quasi-static dense flows (elastic 
and visco-plastic models) are very different, see e.g. 
iGoldenberg and GoldhirschI l)2002|) . However, near 
the fluidization transition, and in dense partially- 
fluidized flows, the differences between these two 
regimes become less obvious. The fluidization of 



sheared granular materials has many features of a 
first-order phase transition. The phenomenologi- 
cal partial fluidization theory in principle can be 
a bridge between the static and dynamic descrip- 
tions. The order parameter related to the local co- 
ordination number appears to be one of the hidden 
fields required for a consistent description of granu- 
lar flows. One important question in this regard is 
the universality of the fluidization transition in dif- 
ferent granular systems and geometries. On the op- 
posite side of the fluidization transition, the static 
state of the granular matter can be described by the 
order parameter related to the percentage of static 
contacts with fully activated dry friction (critical 
contacts) ijStaron et all l2002() . It was shown that 
once these contacts form a percolation cluster, the 
granular pack slips and fluidization occurs. It is of 
obvious interest to relate this "static" order param- 
eter and the "dynamics" order parameter discussed 
above. We see one of the main future challenges in 
the systematic derivation of the continuum theory 
valid both for flowing and static granular matter. 

• Statistical mechanics of dense granular systems. 
Glearly, discrete grain structure plays a major role 
in the dynamics and inherent stochasticity of gran- 
ular response. The number of particles in a typi- 
cal granular assembly is large (10^ or more) but it 
is much smaller than the Avogadro number. Tra- 
ditional tools of statistical physics do not apply 
to dense granular systems since grains do not ex- 
hibit thermal Brownian motion. One of the alter- 
native ways of describ ing statistics of granu l ar me- 
dia was suggested bv lEdwards and GrinevI l)l998|) 
in which they proposed that volume rather than 
energy serves as the extensive variable in a static 
granular system, so that the role of temperature 
is played by the compactivity which is the deriva- 
tive of the volume with res pect to the usual en- 
tropy. Recent expe r iment s ijMakse and Kurchanl 
I2OO2I: ISchroter etdl 120051) aim to test this theory 
experimentally. Gonnecting Edwards theory with 
granular hydrodynamics will be an interesting chal- 
lenge for future studies. 

• Realistic simulations of three-dimensional granu- 
lar flows. Even the most advanced simulat i ons of 
granular flows in t hree dimensions l)Silbertl l2005t 
ISilbert et all l2003|) are limited to relatively small 
samples (e.g. 100 x 40 x 40 particles box) and are 
very time consuming. The granular problems are 
inherently very stiff: while the collisions between 
particles are very short (O(10~'*sec), the collective 
processes of interest may take many seconds or min- 
utes. As a result, to the time step limitations a 
simulation of realistic hard particles is not feasi- 
ble: the "simulations" particles have elastic mod- 
uli several orders of magnitude smaller than sand 
or glass. The particle softness may introduce un- 
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physical artifacts in the overall picture of the mo- 
tion. Different approaches to handling this prob- 
lem will be necessary to advance the state of the 
art in simulations. New opportunity can be offered 
by the equation-free simulation method proposed 
bv iKevrekidis et al\ l|2004l) . Another area of sim- 
ulations which needs further refinements is an ac- 
curate account of dry friction. In the absence of 
a better soluti on current methods (see for review 
ijLudind . l2004j) ') employ various approximate tech- 
niques to simulation dry friction, and accuracy of 
these methods can be questionable. 



• Complex interactions. Understanding of dynam- 
ics of granular systems with complex interactions is 
certainly an intriguing and rapidly developing field. 
While interaction of grains with intersticial fluid is 
a traditional part of engineering research, effects of 
particle anisotropy, long-range electromagnetic in- 
teractions mediating collisions, adhesion, agglom- 
eration and many others constitute a formidable 
challenge for theorists and a fertile field of future 
research. 



• Granular physics on a nano-scale. There is a 
persistent trend in the industry such as powder 
metallurgy, pharmaceutical and various chemical 
technologies towards operating with smaller and 
smaller particles. Moreover, it was recognized re- 
cently that micro- and nano-particles can be use- 
ful for fabrication of desired ordered structures 
and templates for a broad range of nanotech- 
nological applications through self-assembly pro- 
cesses. Self-assembly, the spontaneous organiza- 
tion of materials into complex architectures, con- 
stitutes one of the greatest hopes of realizing the 
challenge to create ever smaller nanostructures. 
It is a particulary attractive alternative to tra- 
ditional approaches such as lithography and elec- 
tron beam writing. Reduction of the particle size 
to micro- and nano -scales shifts the balance be- 
tween forces controlling particle interaction because 
the dominant interactions depend on the parti- 
cle size. While for macroscopic grains the dy- 
namics are governed mostly by the gravity, coUi- 
sional and frictional forces, for micro- and nano- 
particles the dominant interactions include long- 
range electromagnetic forces, short- range van der 
Waals interactions, etc. Nevertheless, some con- 
cepts and ideas developed in the "traditional" 
granular physics were successfully applied to un- 
derstand dynamic s elf-assembly of microparticles 
llSapozhnikov et aZ.l l2003aL 1200 4^ and even biolog- 
ical microtubules ijAranson an d Tsimrine, 2005). 
We expect to see more and more efforts in this di- 
rection. 
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